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1.  Introduction 


As  the  use  of  composite  materials  expands,  the  demand  for  an  accurate  model  of  the  response 
and  failure  of  these  materials  has  increased.  The  adoption  of  the  finite  element  (FE)  method  as  a 
tool  to  analyze  the  response  of  complex  structures  under  load  has  greatly  improved  the 
understanding  and  design  of  these  structures.  Linear,  elastic  materials  that  are  used  in  many 
applications  are  easily  handled  by  any  well-developed  FE  software,  but  the  unique  nature  of 
laminated  composite  materials  have  not  been  fully  addressed  in  commercially  available  software. 
Composite  materials  are  unique  due  to  their  ply-level,  anisotropic,  nonlinear  stress-strain  response. 
It  is  possible  through  the  load-displacement  history  of  a  composite  structure  for  a  ply  or  plies  to  fail 
while  the  overall  structure  does  not  fail.  Accounting  for  first-ply  failure  and  any  subsequent 
failures,  know  as  progressive  failure  analysis,  is  therefore  important  in  accurately  simulating  the 
response  of  composites.  The  incorporation  of  progressive  failure  for  composite  materials  in 
commercial  software  has  been  limited  ( 1 ).  In  this  section,  a  review  of  the  state-of-the-art 
composite  modeling  features  offered  by  the  most  widely  available  FE  software  and  specialized 
composite-oriented  supplemental  software  is  presented.  A  literature  review  on  topics  related  to 
the  progressive  damage  and  failure  analysis  of  composite  structures  is  also  discussed. 

1.1  Problem  Statement 

Large  scale,  thick-section  composite  structures  cannot  be  modeled  efficiently  using  a  discrete 
ply-by-ply  approach.  For  structures  with  hundreds  of  layers,  the  number  of  elements  through  the 
thickness  would  be  too  large  to  be  computationally  efficient.  Creating  the  model  and 
interpreting  the  ply-by-ply  stresses  and  strains  would  be  very  time  consuming  and  prone  to 
errors.  To  address  the  issues  of  ply-by-ply  analysis,  a  multiscale  approach  is  taken  to  account  for 
the  ply-level  response  and  failure  of  the  composite  structure  while  not  explicitly  modeling  the 
plies.  To  achieve  this,  multilayered  laminates  are  treated  as  a  homogenous,  or  smeared,  material 
with  equivalent  material  properties.  Ply-level  material  nonlinearity  and  progressive  failure 
analysis  using  ply-based  failure  criteria  are  incorporated  into  the  multiscale  approach  by 
dehomogenizing,  or  unsmearing,  the  laminate. 

An  investigation  of  three  commercially  available  FE  programs,  state-of-the-art  composite 
modeling  software  products,  and  recent  publications  is  presented  in  this  section.  All  these  works 
contain  either  anisotropic  material  nonlinearity,  or  progressive  failure  analysis,  or  the  multiscale 
“smearing-unsmearing”  approach,  but  it  is  shown  that  none  contain  all  three.  The  combination 
of  these  three  aspects  is  what  makes  the  modeling  approach  presented  in  this  report  unique.  The 
objective  of  this  research  is  to  develop  a  methodology  that  uses  this  unique  multiscale  approach 
in  conjunction  with  FE  analysis  to  design  and  analyze  thick-section  composite  structures. 
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1.2  Benchmarking 

1.2.1  State-of-the-art  Commercial  Composite  Models 

The  most  advanced  composite  modeling  capabilities  of  current,  commercially  available,  FE 
software  are  investigated  in  this  section.  ABAQUS,  ANSYS,  and  MSC.Nastran  are  three  major 
FE  programs  that  have  composite  modeling  capabilities.  All  three  programs  have  specialized 
features  that  enable  the  user  to  construct  and  analyze  composite  materials.  For  the  purpose  of 
this  research,  the  features  investigated  in  these  programs  include  the  ability  to  have  ply-level 
nonlinearity  in  the  principal  and  shear  directions,  progressive  ply  failure  capabilities,  support  for 
a  wide  variety  of  failure  criteria,  the  ability  to  efficiently  model  the  layup  by  homogenizing 
properties  to  get  an  equivalent  response,  and  an  effective  way  of  visualizing  this  data. 

In  recent  software  updates,  ABAQUS  has  incorporated  the  ability  to  model  composite  layups 
through  the  use  of  conventional  or  continuum  shell  layup  sections  for  shell  elements  and  solid 
layup  sections  for  solid  elements  (2).  These  sections  are  able  to  analyze  composite  materials  by 
modeling  each  ply  discretely.  The  layup  can  be  comprised  of  isotropic  nonlinear  materials,  but 
nonlinear  anisotropic  materials  are  not  supported  (3).  There  is  support  for  progressive  failure 
and  damage  evolution  laws  in  ABAQUS,  but  it  is  restricted  to  linear,  elastic  materials,  which  can 
be  anisotropic  (4).  Furthermore,  progressive  failure  behavior  is  restricted  to  shell  elements  and 
the  only  damage  initiation  criterion  that  can  be  used  is  the  Hashin  criteria  (4).  ABAQUS  is  able 
to  visualize  the  damage  initiation  criteria  (tensile  and  compressive  damage  index  for  both  fiber 
and  matrix)  and  has  the  ability  to  label  these  values  by  ply.  Although  limited  in  scope  and 
restricted  to  specific  element  types,  ABAQUS  provides  a  few  important  modeling  features  for 
composite  materials. 

ANSYS  uses  both  shell  (SHEFF99,  91,  and  181)  and  solid  (SOFID186,  46,  191,  and  95) 
elements  to  model  composite  layups  (5).  Each  shell  and  solid  element  listed  has  distinct 
advantages  and  disadvantages,  as  well  as  specific  applications  in  which  they  are  used.  Several  of 
these  elements  are  able  to  handle  nonlinear  material  properties  (5),  with  extensive  options  for 
different  types  of  anisotropic  nonlinear  behavior.  The  three  basic  failure  criteria  that  are 
supported  by  default  are  maximum  stress,  maximum  strain,  and  Tsai-Wu  (<5),  with  the  option  of 
expanding  the  criteria  through  user-written  code.  These  failure  criteria  are  used  to  determine 
whether  the  material  would  have  failed  and  are  not  coupled  to  the  overall  structural  response. 
ANSYS  does  contain  capabilities  for  delamination  and  crack  growth  and  propagation,  but 
progressive  failure  for  composites  in  ANSYS  is  not  possible.  Similar  to  ABAQUS,  ply-based 
results  are  available  for  viewing  in  ANSYS,  but  there  is  limited  support  for  critical  mode  and  ply 
display  options. 

The  programs  developed  under  MSC  Software,  including  Patran,  Nastran,  Marc,  and  several 
others,  have  recently  begun  to  add  progressive  failure  capabilities  to  their  existing  composite 
modeling  abilities.  In  Nastran,  a  total  of  seven  failure  criteria,  including  Puck  and  Hashin,  can 
be  used  in  conjunction  with  progressive  failure  (7).  A  variety  of  elements  are  available  for 
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modeling  composites,  including  two-dimensional  and  three-dimensional  (3-D)  layered  shell, 
solid,  and  solid  shell  elements.  Material  options  can  be  used  to  create  anisotropic,  nonlinear 
material  response.  Like  other  FE  software,  there  are  few  options  for  visualizing  and  analyzing 
composite- specific  design  data  that  is  produced  from  progressive  failure  analysis.  The  majority 
of  MSC  Software’s  progressive  failure  features  are  included  in  a  separate  product,  developed  in 
part  by  the  National  Aeronautics  and  Space  Administration  (NASA),  called  GENOA.  This 
product  and  other  notable  composite  design  oriented  products  are  discussed  in  section  1.2.2. 

1.2.2  State-of-the-art  for  Other  Software  Products 

Before  the  most  recent  releases  of  major  finite  element  analysis  (FEA)  programs,  little  support 
and  few  options  existed  for  modeling  composite  materials.  A  demand  for  accurate 
characterization  of  composite  structures  spurred  the  development  of  several  products  that  went 
beyond  the  available  capabilities  of  these  programs.  Developed  either  as  an  add-on  to  software 
such  as  ABAQUS,  ANSYS,  Nastran,  and  other  programs  or  as  a  standalone  pre-  or  post¬ 
processing  product,  features  such  as  progressive  damage,  nonlinear  material  properties,  and 
advanced  composite  design  visualization  were  able  to  be  incorporated  into  the  FEA  design  cycle. 
While  the  built-in  capabilities  of  current  FEA  programs  are  catching  up,  these  products  still  offer 
some  of  the  most  advanced  composite  modeling  features  available. 

As  mentioned  before,  GENOA  is  standalone  pre-  and  post-processing  software  that  is  compatible 
with  Nastran,  MSC  MARC,  LS-DYNA,  ANSYS,  and  ABAQUS  (8).  Offering  a  multitude  of 
features,  GENOA  includes  multiscale  progressive  failure  analysis,  the  ability  to  handle  nonlinear 
stress-strain  response,  and  advanced  visualization  options.  In  GENOA,  the  plies  that  make  up 
the  composite  structures  are  modeled  and  displayed  separately.  In  the  simulation  of  structures 
that  are  hundreds  of  layers  thick,  the  homogenization  of  composite  laminates  can  improve  the 
computational  efficiency  of  the  model.  A  homogenizing  feature  is  something  that  is  missing 
from  this  software. 

Another  product  that  is  specifically  tailored  to  simulate  composite  materials  is  Helius:MCT 
created  by  Firehole  Technologies.  Helius  was  developed  as  a  user  material  subroutine  for 
ABAQUS  and  ANSYS  (9).  Through  an  ABAQUS  plug-in  and  graphical  interface,  a  composite 
material  that  is  selected  from  a  database  is  turned  into  a  user  material  (UMAT).  This  UMAT 
allows  for  coupled,  progressive  failure  of  the  composite  part  based  on  a  multi-continuum  theory. 
Nonlinear  material  response  is  limited  to  the  in-plane  and  out-of-plane  shear  stiffness  (10). 

Much  like  modeling  a  composite  in  ABAQUS,  a  layered  composite  section  is  created  using  the 
UMAT  as  the  constitutive  material.  The  main  drawback  to  this  approach  is  that  each  ply  is 
modeled  separately,  decreasing  the  computational  efficiency  when  compared  to  a  homogenized, 
equivalent  composite  material.  In  Helius,  the  progressive  failure  analysis  provides  visualization 
of  matrix  and  fiber  failure  of  the  structure,  but  does  not  provide  indication  of  whether  it  is  tensile 
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or  compressive  and  there  is  no  shear  failure  support.  Directionally  dependent  material 
nonlinearity  and  a  homogenization  approach  are  aspects  of  composites  modeling  that  are  not 
currently  available  in  Helius:MCT. 

There  are  several  additional  standalone  and  integrated  products  that  have  been  developed  to 
tackle  a  range  of  composite  modeling  problems.  ESAComp,  offering  integration  with 
ABAQUS,  ANSYS,  LS-DYNA,  Nastran,  NISA,  and  I-DEAS,  is  capable  of  evaluating  several 
different  failure  theories  of  nonlinear  composites  and  provides  informative  plots  of  analysis 
results  through  a  graphical  interface  (11).  NEi  Nastran  (12)  has  unique  composite  analysis 
features  such  as  the  visualization  of  critical  mode  and  ply  data,  progressive  ply  failure  analysis, 
support  for  the  latest  failure  criteria,  and  nonlinear  material  models.  NIS A/Composite  (13)  is 
able  to  display  specialized  composite  analysis  parameters  and  model  nonlinear  materials  using 
several  available  3-D  layered  elements.  The  main  drawback  of  all  these  programs  is  the  reliance 
on  individual  ply  modeling  and  subsequent  lack  of  support  for  the  homogenizing  of  properties, 
which  improves  the  computational  efficiency  of  modeling  of  large-scale,  thick  section, 
composite  structures. 

1.2.3  Literature  Review 

A  review  of  recent  publications  relating  to  progressive  damage  and  failure  modeling,  as  well  as 
nonlinear  composite  material  modeling,  was  conducted  to  determine  the  current  state  of  research 
in  these  areas.  Many  papers  deal  with  both  of  these  areas  and  provide  innovative  ways  to  utilize 
them.  The  search  was  focused  on  trying  to  find  other  research  that  contained  a  part  or  all  of  the 
key  features  that  have  been  determined  to  effectively  model  large  laminated  structures.  The  key 
aspects  of  the  proposed  modeling  solution  are  nonlinear  material  responses  that  are  directionally 
dependent,  progressive  ply  failure,  and  the  ability  to  efficiently  model  the  composite  through 
homogenized,  effective  properties. 

Numerous  researchers  have  conducted  progressive  failure  analyses  using  linear,  orthotropic 
material  properties.  Using  a  newly  developed  damage  evolution  criteria,  Lapczyk  and  Hurtado 
(14)  use  the  progressive  damage  features  within  ABAQUS  to  model  open-hole  tension  tests  of 
aluminum/composite  sandwich  structures.  Zhang  et  al.  (15)  create  a  FE  model  of  grid  stiffened 
plates  subjected  to  damage  based  on  a  modified  Hashin  criterion.  Progressive  damage  and 
delamination  are  investigated  on  open-hole  samples  subjected  to  tension,  compression,  and 
interlaminar  shear.  Spottswood  and  Palazotto  (16)  use  simplified  large  displacement  and 
rotation  theory  to  look  at  the  progressive  failure  of  curved  composite  shells  up  to  and  through 
their  snapping  point.  The  effect  of  two  material  degradation  models  on  the  response  of 
composite  bolted  joints  is  investigated  by  Huhne  et  al.  (17).  Both  a  constant  and  continuous 
degradation  approach  are  incorporated  into  an  ABAQUS  subroutine  and  compared  to 
experimental  data.  Xie  and  Biggers  (18)  look  at  the  effect  width-to-hole-diameter  ratio  on  open- 
hole  tensile  response  by  using  a  progressive  damage  model  that  discounts  plies  instantaneously 
after  failure.  This  behavior  is  programmed  into  an  ABAQUS  user  material  and  the  ability  to 
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tailor  the  composite  using  the  results  is  demonstrated.  Buckling  and  progressive  damage  of  an 
open-hole  compression  sample  is  simulated  by  Labeas  et  al.  (19).  Using  ANSYS,  two  different 
failure  theories  and  two  degradation  models  are  compared  with  experimental  results  to  determine 
the  best  combination.  Tserpes  et  al.  (20,  21)  created  a  user  material  for  ANSYS  that  uses  the 
Hashin  criteria  and  the  ply  discount  method  to  model  progressive  failure.  Variations  on  this 
approach  were  used  to  parametrically  study  laminate -bolt  interactions  to  close  agreement. 

Nonlinear  stress-strain  response,  especially  in  the  transverse  and  through-thickness  directions,  is 
an  important  aspect  of  composites.  There  are  a  number  of  papers  that  incorporate  progressive 
failure  with  nonlinear  response.  Schuecker  and  Pettermann  (22)  create  a  continuum  damage 
model  that  they  compare  to  results  from  the  World  Wide  Failure  Exercise  (WWFE).  Hochard 
et  al.  (23)  develop  an  ABAQUS  user  material  using  shear  damage  and  inelastic  strain  evolution 
laws.  Both  static  and  fatigue  loads  are  applied  to  open-hole  samples  and  strain  fields  are 
compared  to  data  obtained  from  digital  image  correlation.  The  progression  of  failure  and  post- 
buckling  response  of  tapered  composites,  especially  those  with  ply  drop-offs,  is  studied  by 
Ganesan  and  Liu  (24).  Using  nonlinear  ply  properties  and  the  maximum  stress  failure  criteria, 
multiple  layups  are  evaluated  against  one  another.  Kilic  and  Haj-Ali  (25)  developed  an 
ABAQUS  user  material  that  utilizes  Ramberg-Osgood  based  nonlinear  behavior  and  Tsai-Wu 
failure  criteria  to  model  variations  of  both  single  bolt  and  open-hole  tension  tests.  Open-hole  and 
double  notched  compression  samples  are  modeled  by  Basu  et  al.  (26)  using  layered  shell  element 
and  a  user  material  in  ABAQUS.  Nonlinear  material  properties  are  simulated  using  Schapery 
theory  and  fiber  rotation  under  axial  compression.  Antoniou  et  al.  (27)  simulated  thick 
cylindrical  composite  shells  under  compressive,  tensile,  and  torsion  loading  in  ANSYS.  A  user 
material  that  had  linear,  elastic  fiber  direction  response  and  nonlinear  transverse  and  shear 
response  was  used  to  create  failure  envelopes  of  first  and  last  ply  failure  for  comparisons  with 
experimental  data.  Huang  (28)  developed  an  ABAQUS  user  shell  section  that  uses  the  Bridging 
model  to  degrade  the  nonlinear  fiber  and  matrix  properties  upon  failure.  The  effect  of  mesh  size, 
element  type,  and  load  application  are  investigated  and  compared  to  experimental  results.  A 
unique  aspect  of  this  failure  analysis  is  that  if  a  ply  associated  with  an  element  fails,  that  ply  is 
discounted  throughout  the  entire  structure. 

While  there  are  several  examples  of  both  linear  and  nonlinear  behavior  coupled  with  progressive 
failure  analysis,  none  of  these  papers  have  the  ability  to  simulate  effective  properties  by 
homogenizing  the  laminate  into  one  material.  An  important  part  of  modeling  is  the 
computational  efficiency  of  the  code,  so  the  designer  can  quickly  analyze  many  design  iterations. 
This  aspect  of  progressive  failure  and  composites  simulation  has  not  been  addressed  in  previous 
work  and  is  important  when  creating  an  effective  design  tool. 

1.2.4  Selection  of  FEA  Software 

An  investigation  of  the  phrases  composite,  nonlinear,  nonlinear  composite,  composite  progressive 
failure,  and  nonlinear  composite  progressive  failure  was  conducted  using  publications  from 
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January  2000  to  January  2010.  The  most  widely  mentioned  codes  in  publications  containing  the 
key  phrases  were  ABAQUS,  ANSYS,  and,  to  a  lesser  extent,  MSC.Nastran.  Table  1  shows  the 
percentage  of  papers  with  the  phrases  that  also  mention  ABAQUS,  ANSYS,  or  Nastran. 
Throughout  all  of  these  papers,  the  most  widely  used  program  is  ABAQUS,  ranging  from  1%  of 
all  papers  with  the  keyword  composite  to  17.6%  of  all  papers  with  the  keywords  nonlinear 
composite  progressive  failure. 

Table  1.  Percentage  of  keyword  papers  containing  FEA  software. 


Keywords 

Abaqus 

(%) 

ANSYS 

(%) 

Nastran 

(%) 

Composite 

1.1 

0.8 

0.2 

Nonlinear 

1.3 

0.8 

0.1 

Nonlinear  composite 

4.7 

3.0 

0.5 

Composite  progressive  failure 

5.6 

2.5 

0.4 

Nonlinear  composite  progressive  failure 

17.6 

6.8 

1.3 

Due  to  its  wide  use  in  modeling  the  progressive  damage  and  failure  of  nonlinear  composite 
materials,  ABAQUS  was  selected  as  the  software  in  which  the  homogenization  approach  and 
progressive  ply  failure  technique  will  be  implemented. 


2.  Solid  Mechanics  Approach 


The  solid  mechanics  based  “LAM”  codes  were  developed  for  the  design  and  failure  analysis  of 
thick-section  composite  structures.  Both  linear  and  nonlinear,  these  codes  use  a  multiscale 
“smearing-unsmearing”  approach  and  progressive  ply  failure  to  analyze  multilayered  composite 
laminates.  A  background  on  the  development  of  these  programs  and  discussion  on  the  key 
assumptions  and  procedures,  particularly  for  the  nonlinear  codes,  are  given  in  this  section. 

2.1  Smearing-Unsmearing  Approach 

The  analysis  of  composite  structures  using  FE  techniques  is  complicated  by  the  nonlinear, 
anisotropic  ply-level  response  of  these  materials.  Furthermore,  the  failure  analysis  of  a  ply  or 
lamina  depends  upon  the  stress  and  strain  state  within  individual  plies.  It  is  therefore  important 
to  track  the  individual  ply  stresses  and  strains  in  order  to  accurately  model  the  nonlinear  material 
response  and  to  be  able  to  apply  failure  predictions. 

One  way  to  accurately  calculate  the  ply-level  stresses  and  strains  would  be  to  model  each  ply 
discretely.  Using  several  elements  through  the  thickness  of  the  ply,  a  mesh  containing  all  the 
plies  of  the  structure  could  be  created.  While  this  approach  would  most  accurately  capture  the 
ply-level  response,  there  would  be  several  disadvantages.  A  model  that  is  both  large  in  scale  and 
in  number  of  plies  would  prove  very  time  consuming  to  create,  mn,  and  analyze.  In  addition,  the 


6 


vast  amount  of  data  produced  by  a  large,  complex  model  would  be  difficult  to  interpret.  There 
would  also  be  a  greater  chance  for  errors  when  creating  the  model  when  compared  to  a 
homogenized  structure. 

A  way  to  address  the  computational  and  logistical  disadvantages  of  discretely  modeling  each  ply 
is  by  homogenizing  the  laminate.  Complications  arise  when  the  application  of  ply-based  failure 
criteria  are  considered  for  a  homogeneous  material.  Ply-level  stresses  and  strains  must  be 
extracted  from  structural  stress  and  strain  profile  in  order  to  conduct  failure  analysis.  To  deal 
with  these  issues,  Bogetti  et  al.  implemented  the  well-established  “smearing/unsmearing” 
approach  (29-31).  Equivalent  laminate  properties  are  generated  by  “smearing”  the  ply-level 
nonlinear,  anisotropic  material  behavior  of  the  layup  and  ply-level  stress  state  needed  for  failure 
analysis  is  generated  by  “unsmearing”  the  elemental  stresses  and  strains. 

2.2  Background  of  the  “LAM”  Codes 

The  programs  LAM3D,  LAM3DNLP,  LAMPAT,  and  LAMPATNL  are  the  four  codes  in  the 
“LAM”  series.  The  programs  use  the  “smearing-unsmearing”  approach  and  progressive  ply 
failure  analysis  to  model  multilayered  composite  laminates  on  two  scales.  LAM3D  and  its 
nonlinear  version  LAM3DNLP  are  codes  that  use  analytical  models  to  simulate  the  material 
point  response  of  laminates.  LAMPAT  and  LAMPATNL  scale  up  the  LAM3D  and 
LAM3DNLP  approaches  for  use  in  FE  analysis  to  evaluate  composite  structures.  A  discussion 
of  each  code  is  provided  in  this  section. 

LAM3D  was  created  by  Bogetti  et  al.  (29)  and  uses  orthotropic,  linear  ply  properties  and 
progressive  ply  failure  to  simulate  the  material  point  response  of  a  laminate.  The  code 
determines  the  3-D  effective  properties,  ultimate  strength  predictions  under  mechanical  and 
thermal  loading,  and  effective  stress-strain  response  due  to  progressive  ply  failure  under 
mechanical  loading  of  thick-section  composites.  Using  a  specified  mechanical  load,  LAM3D 
simulates  first  ply  to  last  ply  failure,  calculates  effective  laminate  properties  and  ply-level 
stresses  and  strains  after  each  failure,  and  determines  the  safety  factor,  critical  mode,  and  critical 
ply  of  the  laminate. 

The  effective  property  and  strength  predictions  of  the  LAM3D  code  for  many  composite 
laminates  were  validated  against  theoretical  predictions  and  experimental  results.  While 
LAM3D  provided  valuable  insight  into  laminate  response  due  to  progressive  failure,  the  code  did 
not  support  nonlinear  ply  properties.  Nonlinear  material  response  is  more  accurate  than  linear 
for  the  transverse  and  shear  directions  of  composite  materials.  In  order  to  address  this  issue, 
Bogetti  et  al.  (30)  created  the  nonlinear,  progressive  failure  code  LAM3DNLP.  In  this  code, 
nonlinear  constitutive  relationships  for  all  six  principal  directions  were  modeled  using  Ramberg- 
Osgood  parameters.  Similar  to  LAM3D,  the  “smearing-unsmearing”  approach  and  progressive 
ply  failure  methodology  are  used  by  LAM3DNLP  to  simulate  the  material  point  response  of  a 
laminate.  Both  LAM3D  and  LAM3DNLP  were  independent  analytical  programs  developed 
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using  FORTRAN  and  were  not  coupled  to  FE  analysis.  The  desire  to  use  these  material  point 
models  in  the  design  of  large  composite  structures  led  to  the  development  of  LAMPAT  and 
LAMPATNL. 

The  LAM3DNLP  analytical  code  was  ranked  third  out  of  19  participants  in  the  first  World  Wide 
Failure  Exercise  (WWFE-I)  (32).  This  exercise  was  an  assessment  of  the  state  of  the  art  of 
composite  failure  predictions.  LAM3DNLP  was  used  to  predict  progressive  ply  failure  in 
composite  laminates  for  14  loading  cases.  These  failure  predictions  were  compared  to 
experimental  results  and  rated  against  the  failure  predictions  of  the  other  participants  (32-34). 
The  14  loading  cases  provided  an  experimental  validation  of  the  LAM3DNLP  failure 
predictions.  The  first  exercise  focused  on  failure  response  of  laminates  subjected  to  in-plane 
loading.  The  second  exercise  (WWFE-II)  is  currently  in  progress  (35).  This  exercise  focuses  on 
out-of-plane  failure  predictions  and,  similar  to  the  first  exercise,  will  provide  additional 
validation  for  LAM3DNLP. 

LAMPAT  was  created  by  Bogetti  et  al.  (36)  as  a  way  of  using  the  “smearing-unsmearing” 
approach  and  progressive  failure  analysis  of  LAM3D  with  FE  software  packages.  LAMPAT 
was  designed  as  a  pre-  and  post-processor  that  created  effective  composite  laminate  properties 
and  evaluated  the  failure  of  laminates  based  on  element  stress  and  strain.  Originally  developed 
to  interface  with  PATRAN,  it  was  later  expanded  to  accommodate  ANSYS,  ABAQUS, 
DYNA3D,  and  NIKE3D.  As  a  pre -processor,  LAMPAT  uses  the  orthotropic,  linear  lamina 
properties  and  the  architecture  of  the  laminate  to  create  a  set  of  “smeared”  material  properties 
that  represent  the  effective  response  of  the  laminate.  Once  the  composite  structure  with  effective 
properties  is  evaluated,  the  postprocessor  takes  global  stresses  and  strains  and  creates 
“unsmeared”  ply  stresses  and  strains.  The  stress  and  strain  allowables  of  the  lamina  are  then 
used  to  analyze  failures  within  the  composite  to  provide  safety  factor,  critical  mode,  and  critical 
ply  data.  This  process  is  shown  in  figure  1. 

Implementing  LAM3D  into  FE  software  using  LAMPAT  demonstrated  the  power  of  this 
approach  when  designing  large-scale,  thick-section  composite  structures.  To  more  accurately 
capture  the  response  of  composite  materials,  the  material  nonlinearity  features  of  LAM3DNLP 
were  needed.  In  LAMPAT,  failure  analysis  occurred  after  the  stress  and  strain  profile  of  the 
structure  was  complete.  This  did  not  allow  for  the  redistribution  of  load  through  the  structure 
after  ply  failure  had  occurred.  In  order  to  address  this  issue,  Powers  et  al.  (36)  developed  the 
nonlinear,  progressive  failure  code  LAMPATNL.  While  LAMPAT  was  developed  as  a  pre-  and 
post-processing  tool,  the  nonlinear  material  behavior  and  progressive  failure  analysis  of 
LAMPATNL  required  closer  integration  to  the  FE  package.  This  was  accomplished  by  creating 
a  coupled  analysis  using  a  user-defined  material  in  ABAQUS.  Limited  validation  and 
demonstration  of  the  LAMPATNL  code  was  provided  by  Powers  in  (36).  A  full  validation  and 
development  of  a  design  methodology  using  LAMPATNL  is  documented  in  sections  3  and  4. 
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Figure  1 .  LAMP  AT  process  methodology. 


2.3  Laminated  Media  Methodology 

The  3-D  laminated  media  model  that  is  the  basis  of  all  the  “LAM”  codes  was  created  from  an 
analytical  model  developed  by  Chou  et  al.  (38).  The  model  is  used  to  determine  the 
homogenized,  or  “smeared,”  engineering  constants  of  a  multilayered  laminate  (figure  2). 
Significant  features  of  the  theory  are  discussed  by  Bogetti  et  al.  in  (29-31)  and  an  overview  is 
given  in  this  section. 
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Equation  1  represents  the  effective  stress/strain  constitutive  relationship  for  the  laminate: 
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The  barred  notation  is  used  to  identify  that  it  is  in  the  coordinate  system  of  the  laminate  and  the 
star  (*)  signifies  that  it  is  an  average  value.  The  relationship  between  the  global  coordinates 
(X,  Y,  Z),  laminate  coordinates  (  X  5  Y ;  Z  );  and  ply  coordinates  (1,  2,  3)  is  shown  in  figure  3. 


Figure  3.  Relation  of  global  coordinates  to  laminate  and  ply  coordinates. 

Chou  et  al.  (37)  derived  a  homogeneous  representation  for  the  laminate.  The  coefficients  of  the 
laminate  stiffness  matrix,  Cj,  are  given  in  equations  5-7. 
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Where  the  k  term  refers  to  the  k,b  ply  of  the  laminate,  4  is  the  ratio  of  the  original  thickness  of 
the  kth  ply  over  the  original  thickness  of  the  entire  laminate,  and 
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The  assumption  is  made  that  the  applied  mechanical  loading  acting  on  the  laminate  [cr  J  is 
known,  not  varying  through  the  thickness,  and  represents  the  “average”  or  “effective”  stress. 
The  associated  “effective”  or  “smeared”  laminate  strains  [/7  ]  can  be  obtained  directly  from  the 
inversion  of  equation  1 . 


In  determining  the  individual  ply-level  stresses  and  strains,  two  major  assumptions  are  used. 
These  assumptions  are  that  the  in -plane  ply  strains  are  equal  to  the  effective  strains  of  the 
laminate  and  that  the  out-of-plane  ply  stresses  are  uniform  through  the  thickness  and  equal  to 
the  effective  stresses  in  the  laminate.  The  assumptions  are  shown  in  equations  9  and  10, 
respectively. 

£■=£■  for(i=  1,2,  6andk=  1,2,  ...,N)  (9) 


<jk  =  a*  for  (i  =  3,  4,  5  and  k  =  1,2,  . . . ,  N) 


(10) 


The  expression  derived  by  Sun  and  Liao  (38)  is  used  to  determine  the  remaining  strain 
components  and  is  given  in  equation  11. 
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The  remaining  stress  components  are  easily  calculated  using  the  relation  in  equation  12. 
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for  (k  =  1,2,  ...,  N) 


(12) 


2.4  Nonlinear  Constitutive  Relationships 

Ply-level  material  nonlinearity  is  represented  by  using  the  Ramberg-Osgood  equation.  Three- 
dimensional  nonlinearity  is  accommodated  in  all  principal  material  directions.  The  longitudinal, 
transverse,  through  thickness,  out-of-plane  shear,  and  in-plane  shear  directions  can  all  have 
distinct  nonlinear  behavior.  The  Ramberg-Osgood  equation,  shown  below,  defines  stress  as  a 
function  of  strain  and  three  parameters. 


Eqg 
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(  f  -n.  VV 

1+  E°S 
l  v  ^0  J  J 

(13) 

where  E0  is  the  initial  modulus,  a(l  is  the  stress  asymptote,  and  n  is  the  shape  factor  of  the 
nonlinear  stress-strain  curve.  The  effect  that  these  parameters  have  on  the  stress-strain  response 
is  shown  in  figure  4. 
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Figure  4.  Nonlinear  stress-strain  response. 


For  computational  considerations,  it  is  desired  to  define  the  instantaneous  or  tangent  lamina 
stiffness  as  a  continuous  function  of  strain.  Taking  the  derivative  of  equation  13  with  respect  to 
strain,  the  following  expression  is  obtained: 
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where  Et  is  the  instantaneous  or  tangent  lamina  stiffness  modulus  expressed  explicitly  in  terms  of 
strain  and  the  three  Ramberg-Osgood  parameters. 

2.5  Incremental  Approach  to  Predicting  Nonlinear  Laminate  Response 

The  analytical  code  LAM3DNLP  models  the  nonlinear  laminate  response  as  a  set  of  piecewise 
linear  increments.  It  predicts  laminate  stress-strain  response  using  a  specified  load  increment  up 
until  failure.  The  incremental  approach  starts  with  the  initial  stiffness  as  determined  by  the 
smeared  properties  of  the  laminate.  At  the  end  of  each  load  increment,  the  stress  state  from  the 
laminate  is  unsmeared  for  ply-level  stresses  and  strains,  failure  analysis  occurs,  the  nonlinear 
constitutive  relationships  are  updated  to  determine  the  stiffness  in  each  principal  ply  direction, 
and  effective  properties  used  for  the  application  of  the  next  load  increment  are  obtained  by 
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smearing  the  laminate.  Figure  5  provides  a  representation  of  the  incremental  loading  strategy  for 
an  arbitrary  laminate. 


Figure  5.  Incremental  laminate  loading  methodology  (30). 


ponding  to  the  end  of  the  n'h  stress  increment,  the  strain  and  str 
[<£*"  ,  rx  "  ]).  From  this  point,  the  objective  is  to  determine  the 

L  \ _ /  —  1  f—  *H+ll\  rrTi_  _  _rr _ .* _ 1 _ • _ .  _ _ . 


Assume  that  at  point  (a),  corresponding  to  the  end  of  the  n  stress  increment,  the  strain  and  stress 
state  of  the  laminate  is  known  ( | 
strain  and  stress  state  at  point 
the  end  of  stress  increment  n, 
media  model  constitutive  relation,  equation 
corresponding  increment  in  laminate  strain, 
equation  1: 


fb)  or  ( \s  *n ' 1  ],  [cr  n+1  ]).  The  effective  laminate  stiffness  matrix  at 
\C*n  ,  is  computed  from  an  incremental  form  of  the  laminated 

L  1  J  r  — *  i 

.  With  the  increment  in  load  defined,  [Acr;.  "  J,  the 
As  "  ],  is  calculated  from  an  inverse  form  of 


(15) 


Individual  ply  stress  and  strain  increments  are  calculated  according  to  the  equations  presented 
previously.  A  cumulative  summation  is  maintained  to  track  the  total  stress-and-strain  levels  in 
each  ply  of  the  laminate.  The  tangent  modulus  values  for  each  ply  and  material  direction  are 
calculated  according  to  equation  14  and  used  in  the  determination  of  the  laminate  stiffness  matrix 
for  the  next  laminate  stress  increment  calculation.  The  entire  nonlinear  response  for  the  laminate 
is  obtained  by  the  cumulative  sum  of  all  stress  and  strain  increments  throughout  the  entire  stress 
loading  history. 
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Figure  6  illustrates  the  importance  of  the  incremental  approach  when  modeling  nonlinear 
material  behavior.  This  figure  shows  the  modeled  lamina  stress-strain  response  with  three 
different  load  increments  applied.  The  solid  line  represents  the  exact  response  that  is  the 
analytical  solution  of  the  Ramberg-Osgood  equation.  The  larger  load  increment  results  in  a 
response  that  is  stiffer  than  the  desired  nonlinear  response.  As  the  increment  size  decreases,  the 
model  converges  towards  the  exact  response.  It  is  therefore  important  to  select  a  sufficiently 
small  load  increment  when  simulating  material  nonlinearity. 


Figure  6.  Effect  of  increment  size  on  nonlinear  stress-strain  response. 

2.6  Progressive  Ply  Failure  Methodology 

Along  with  ply-level  nonlinearity,  progressive  ply  failure  will  also  influence  the  effective 
laminate  response  to  load.  Using  the  stresses  and  strains  of  individual  plies,  failure  predictions 
and  post-failure  behavior  are  modeled  for  the  laminate.  Failure  criteria  that  have  been 
incorporated  into  the  “LAM”  codes  include  Maximum  Strain,  Maximum  Stress,  Hydrostatic 
Pressure  Adjusted  Maximum  Stress,  Tsai-Wu,  Christensen,  Feng,  Hashin,  and  Von  Mises.  The 
maximum  strain  failure  criterion  is  commonly  used  for  analysis  using  the  “LAM”  codes  (34,  36). 
This  criterion  compares  the  current  strains  of  each  ply  (si,  82,  S3,  84,  85,  86)  to  the  maximum 
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strain  allowables  (Y1T,  Y1C,  Y2T,  Y2C,  Y3T,  Y3C,  Y23,  Y13,  Y12)  in  each  principal  direction 
and  assumes  a  ply  fails  if  these  strain  allowables  are  exceeded.  The  relations  of  the  strains  to 


allowables  are  documented  in  equations  16-24. 

If  £]>0  and  if  8i>Y IT,  then  the  failure  mode  is  fiber  tension  (16) 

If  £i<0  and  if  l£i l>Y  1C,  then  the  failure  mode  is  fiber  compression  (17) 

If  £2>0  and  if  £2>Y2T,  then  the  failure  mode  is  matrix  tension  (18) 

If  £2<0  and  if  l£2l>Y2C,  then  the  failure  mode  is  matrix  compression  (19) 

If  £3>0  and  if  £3>Y3T,  then  the  failure  mode  is  matrix  tension  (20) 

If  £3<0  and  if  l£3l>Y3C,  then  the  failure  mode  is  matrix  compression  (21) 

If  l£4l>Y23,  then  the  failure  mode  is  interlaminar  shear  (22) 

If  l£sl>Y  13,  then  the  failure  mode  is  interlaminar  shear  (23) 

If  l£6l>Y12,  then  the  failure  mode  is  in-plane  shear  (24) 


The  ply  discount  method  is  used  by  the  “LAM”  codes  to  reduce  the  material  properties  of  failed 
plies.  In  this  method,  the  material  stiffness  in  the  direction  of  failure  is  discounted  and  the 
Poisson’s  ratios  associated  with  that  direction  are  adjusted.  For  example,  if  a  ply  fails  in  the  2 
direction,  the  elastic  modulus  of  that  ply  in  the  2  direction  is  set  to  1%  of  the  current  value  and 
this  direction  is  decoupled  by  setting  V|2  and  V23  to  zero.  Decoupling  the  ply  by  zeroing  the 
Poisson’s  ratios  prevents  unrealistic  ply  strains  caused  by  Poisson’s  effect  in  the  failed  ply. 
Table  2  lists  which  properties  are  modified  for  the  six  different  modes  of  failure. 


Table  2.  Failure  mode  and  property  modification. 


Failure  Mode 

Property  Modification 

1 

E[  =  0.01*Eb  V12  =  v13  =  0 

2 

E2  =  0.01  *E2,  V12  =  v23  =  0 

3 

E3  =  0.01  *E3,  vn  =  v23  =  0 

12 

G12  =  0.01*G12 

13 

g13  =  o.oi*g13 

23 

g23  =  0.01  *g23 

Bogetti  has  described  the  progressive  failure  procedure  of  the  analytical  model  (30).  As  the 
laminate  is  loaded  and  laminate  strains  develop,  the  individual  ply  strains  are  monitored.  When 
ply  failure  is  predicted  in  any  ply,  according  to  the  maximum  strain  failure  criteria,  the 
incremental  loading  to  that  point  is  stopped  and  the  entire  laminate  stress  versus  strain  response 
is  recorded.  The  modulus  associated  with  the  particular  mode  of  failure  in  the  failed  ply  is  then 
reduced  according  to  property  modifications  in  table  2  and  the  incremental  loading  strategy  is 
repeated  from  the  beginning  (all  stresses  and  strains  are  set  to  zero).  The  loading  procedure  is 
continued  until  the  next  failure  in  a  ply  is  detected.  The  corresponding  modulus  value  is  again 
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discounted,  the  laminate  response  is  recorded,  and  the  procedure  is  repeated.  This  progressive 
ply  failure  response  is  repeated  until  final  failure  is  determined,  which  is  assumed  when  the 
laminate  loses  sufficient  stiffness  such  that  it  cannot  carry  any  load  without  undergoing  an 
arbitrarily  excessive  amount  of  deformation  (>10%  strain). 

This  progressive  failure  procedure  is  effective  in  capturing  the  nonlinear  stress-strain  response  of 
the  laminate,  but  modifications  to  this  approach  are  needed  for  successful  implementation  into 
FEA  software.  These  modifications  are  discussed  in  section  3. 


3.  Implementation 


LAM3DNLP  was  implemented  into  FE  software  by  Powers  et  al.  with  the  creation  of 
LAMPATNL  (4).  LAMPATNL  is  a  user  material  subroutine  (UMAT)  used  to  incorporate  the 
“smearing-unsmearing”  approach,  nonlinear  anisotropy,  and  progressive  failure  analysis  into 
ABAQUS.  The  subroutine  UMAT  is  used  to  define  the  mechanical  constitutive  behavior  of  a 
material  and  is  called  at  all  material  calculation  points  of  the  element.  LAMPATNL  treats  each 
material  calculation  point  within  the  FE  model  as  though  it  was  a  laminate  in  LAM3DNLP. 

Several  implications  arise  when  converting  the  analytical  code  into  a  UMAT  subroutine.  The 
analytical  code  was  designed  to  simulate  a  laminate  as  a  material  point  under  a  prescribed  load 
increment.  While  it  is  relatively  straightforward  to  keep  track  of  the  nonlinear  response  and 
failure  history  of  a  single  point,  mapping  this  procedure  for  thousands  of  elements  and  material 
points  can  be  complicated.  In  addition,  there  are  certain  constraints  imposed  by  the  FE  method 
to  insure  stability.  The  constraints  are  based  on  thermodynamic  stability  conditions.  Changes  to 
the  material  response  and  the  key  assumptions  of  the  laminated  media  methodology  are  required. 
Powers  (36)  briefly  describes  these  implications  and  changes,  but  an  investigation  into  these 
issues  is  discussed  below. 

3.1  Theoretical  Consideration 

3.1.1  Ply-level  Stress  and  Strain  Calculations 

The  two  major  assumptions  of  the  theory  used  in  the  “LAM”  codes  are  stated  in  section  2.  The 
first  assumption  is  that  the  in-plane  ply  strains  are  equal  to  the  effective  strains  of  the  laminate 
and  is  shown  in  equation  25.  The  second  assumption  is  that  the  out-of-plane  stresses  are  uniform 
and  equal  to  the  effective  stresses  in  the  laminate  and  is  expressed  in  equation  26.  In  these 
equations,  the  k  term  refers  to  the  klh  ply  of  the  laminate. 
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£■=£*  for  (i  =  1,2,  6andk=  1,2,  ...,N)  (25) 

crk  =  a*  for  (i  =  3,  4,  5  and  k  =  1,2,  . . .,  N)  (26) 

All  the  remaining  ply  strains  (3,  4,  and  5)  and  ply  stresses  (1,2,  and  6)  are  determined  by  the 
expressions  given  by  Sun  and  Liao  (3<8),  shown  in  equations  27  and  28,  respectively. 
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for  (k  =  1,2,  ...,  N) 


(28) 


The  ply  stresses  and  strains  are  used  to  calculate  the  tangential  stiffness  properties  for  each  ply. 
These  ply-level  properties  are  “smeared”  to  determine  the  effective  stiffness  of  the  laminate  for 
that  increment.  Ply  stresses  and  strains  are  also  needed  to  apply  ply-based  failure  criteria  for 
progressive  failure  in  the  structure. 

The  procedure  for  calculating  the  ply-level  stresses  and  strains  for  the  nonlinear  analytical  code 
is  outlined  in  figure  7.  In  LAM3DNLP,  the  3-D  load  increment  for  the  laminate,  in  the  form 
( cr  r ,  <x ,  cr. ,  <x  ,  (txz  ,  cr  ),  is  input  into  the  program.  After  initializing  laminate  strain,  ply 
strains,  and  ply  stresses,  the  program  calculates  the  ply-level  tangential  stiffness  properties, 
laminate  stiffness  matrix,  and  ply-level  stiffness  matrices  based  upon  the  current  ply  strains.  The 
laminate  stiffness  matrix  and  laminate  load  increment  are  then  used  to  calculate  the  laminate 
strain  increment.  Using  the  assumption  in  equation  25,  the  in-plane  ply  strain  increment  is  set 
equal  to  the  in-plane  laminate  strain  increment.  Similarly  for  the  assumption  in  equation  26,  the 
out-of-plane  ply  stress  increment  is  set  equal  to  the  out-of-plane  laminate  stress  increment.  The 
out-of-plane  ply  strain  increment  is  calculated  using  the  ply-level  stiffness  matrices,  the  in-plane 
laminate  strain  increment,  and  the  out-of-plane  laminate  stress  increment  using  equation  27.  The 
in-plane  ply  stress  increment  is  calculated  similarly  using  equation  28.  The  increments  for  ply 
stress,  ply  strain,  and  laminate  strain  are  used  to  update  their  respective  values  and  the  procedure 
is  continued  for  the  next  load  increment.  This  procedure  demonstrates  the  direct  implementation 
of  the  assumptions  of  Chou’s  theory  when  calculating  the  ply- level  stress  and  strain  components 
in  LAM3DNLP. 
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Known  value  from  input: 


Figure  7.  Ply  stress  and  strain  calculation  procedure  for  LAM3DNLP. 
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Implementing  LAM3DNLP  into  the  LAMPATNL  UMAT  requires  modifications  to  the 
procedure  shown  in  figure  7.  In  ABAQUS,  the  UMAT  is  required  to  output  the  stiffness  matrix 
and  update  the  stress  tensor  at  the  material  calculation  point  at  the  end  of  the  increment.  These 
values  must  be  determined  from  the  current  laminate  stress,  the  current  laminate  strain,  and  the 
laminate  strain  increment.  The  procedure  for  calculating  the  stiffness  matrix  and  stress  tensor 
using  ply-level  stresses  and  strains  in  LAMPATNL  is  outlined  in  figure  8. 

The  major  difference  between  the  LAMPATNL  procedure  and  the  LAM3DNLP  procedure 
occurs  when  determining  the  tangential  stiffness  properties  of  the  plies.  In  LAM3DNLP,  the  ply 
strains  that  are  initially  zero  are  recorded  and  updated  throughout  the  procedure.  These  strains 
are  used  to  calculate  ply-level  material  properties  and  ply  and  laminate  stiffness  matrices.  The 
LAMPATNL  user  material  subroutine  does  not  record  the  ply-level  strain  data  from  increment  to 
increment  due  to  the  “smearing”  approach  in  the  code.  LAMPATNL  must  calculate  the  ply- 
level  material  properties  and  ply  and  laminate  stiffness  matrices  at  the  beginning  of  the 
increment  without  the  ply-level  strains.  To  accomplish  this,  the  ply  strains  are  set  equal  to  the 
laminate  strains  as  shown  in  equation  29. 

£■=£*  for(i=  1,2,  3,4,  5,  6andk=  1,2,  ...,N)  (29) 

This  enables  LAMPATNL  to  calculate  the  ply  and  laminate  stiffness  matrices  at  the  beginning  of 
the  increment,  update  the  stress  tensor  for  the  increment,  and  determine  the  laminate  stiffness 
matrix  at  the  end  of  the  increment.  Ply  stresses  used  for  stress-based  failure  criteria  are 
determined  analytically  using  the  ply  strains  and  Ramberg-Osgood  relationship. 

The  out-of-plane  ply  strains  are  taken  directly  from  the  laminate  strains  in  LAMPATNL  while 
the  out-of-plane  strains  in  LAM3DNLP  are  calculated  using  the  stiffness  matrix,  out-of-plane 
laminate  stresses,  and  in-plane  laminate  strains.  These  changes  do  not  affect  the  calculated 
stress-strain  response  of  laminate  made  of  one  material  through  the  thickness  or  a  laminate 
subjected  to  in-plane  loading.  The  implication  of  this  change  for  hybrid  composites,  with 
different  material  through  the  thickness,  is  a  subject  for  future  work. 
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Figure  8.  Ply  stress  and  strain  calculation  procedure  for  LAMPATNL. 
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3.1.2  Poisson’s  Ratio 


The  response  of  the  composite  laminate  is  modeled  using  the  Ramberg-Osgood  equation  given  in 
equation  13.  This  stress-strain  response  is  nonlinear  and,  when  coupled  with  the  incremental 
approach,  becomes  a  piecewise  linear  response.  At  each  load  step,  the  instantaneous  tangent 
moduli  for  the  extensional  and  shear  directions  are  calculated  using  equation  14.  Over  each  load 
step,  the  material  is  considered  to  be  a  linear,  elastic  material. 


The  numerical  stability  for  calculations  of  a  linear,  elastic,  orthotropic  material  requires  the 
following  five  conditions  to  be  met  (40): 
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At  the  start  of  the  analysis,  these  conditions  are  satisfied,  but  the  incorporation  of  ply-level 
progressive  failure  can  cause  conditions  31-33  to  become  unsatisfied.  To  correct  this  problem, 
the  Poisson’s  ratio  for  each  ply  is  continuously  updated  throughout  the  simulation.  The  current 
moduli  are  used  in  conjunction  with  the  original  minor  Poisson’s  ratios  (V21,  V31,  V32),  which  are 
stored  at  the  beginning  of  the  analysis,  to  calculate  the  updated  Poisson’s  ratios  and  ensure 
stability.  The  reciprocity  relationships  used  for  this  are 
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where  v  is  the  current  Poisson’s  ratio,  v  is  the  initial  Poisson’s  ratio,  and  ET  is  the 
tangential,  or  instantaneous,  elastic  modulus  at  the  load  increment.  These  adjustments  ensure 
that  the  material  stability  of  the  UMAT  is  satisfied  through  the  entire  simulation.  The 
LAMPATNL  work  of  Powers  et  al.  ( 4 )  contains  these  adjustments  but  does  not  specifically 
document  them. 

3.2  Procedural  Considerations 

3.2.1  Progressive  Failure  Adaption  for  FEA 

The  LAM3DNLP  approach  was  developed  to  analyze  the  material  point  response  of  a  laminate. 
Modifications  to  the  post  failure  behavior  of  the  analysis  are  made  when  implementing  the  code 
in  LAMPATNL.  In  LAM3DNLP,  a  specified  load  increment  is  applied  to  the  laminate  up  until 
failure.  Once  failure  occurs,  the  necessary  adjustments  to  the  failed  plies,  such  as  elastic 
modulus  reduction  and  Poisson’s  ratio  adjustment,  are  made  and  the  analysis  of  the  modified 
material  point  restarts  from  the  first  load  increment.  This  process  is  illustrated  in  the  flowchart 
shown  in  figure  9. 


Figure  9.  LAM3DNLP  progressive  failure  procedure. 
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The  procedure  creates  multiple  stress-strain  curves  and  load  “drops”  between  them,  as  shown  in 
figure  10.  This  figure  shows  the  stress-strain  response  and  progressive  failure  for  this  quasi¬ 
isotropic  layup  of  S2  Glass/Epoxy.  First  ply  failure  is  represented  in  the  first  iteration.  Ultimate 
or  last  ply  failure  is  defined  as  the  failure  that  occurs  at  the  highest  stress  state.  In  this  case,  the 
third  iteration  is  determined  to  be  representative  of  last  ply  failure. 


Figure  10.  Stress-strain  curves  of  [0°/90o/+45°]s  S2  glass/epoxy  laminate. 

The  effective  stress-strain  response  of  the  laminate  is  the  combination  of  the  first  three  iterations. 
This  output  is  represented  as  the  circled  LAM3DNLP  stress-strain  response  in  figure  10.  This 
iterative  procedure  is  effective  in  capturing  the  immediate  loss  of  load  carrying  capability  that  is 
seen  in  experimental  stress-strain  curves,  but  the  process  of  restarting  the  analysis  from  the  first 
increment  is  not  practical  for  FE  applications. 

Changing  this  approach  when  implementing  FAMPATNF  was  done  simply  by  continuing  the 
analysis  with  reduced  material  properties  at  the  current  stress-strain  state  after  failure  has 
occurred.  This  change  is  reflected  in  the  FAMPATNF  flowchart  of  figure  11.  Comparing  this 
procedure  to  FAM3DNFP  procedure  in  figure  9,  the  stress  and  strain  state  of  the  laminate  is  no 
longer  reinitialized  and  laminate  is  not  restarted  from  the  first  increment. 
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Figure  1 1 .  LAMPATNL  progressive  failure  procedure. 


The  process  changes  the  predicted  stress-strain  response  to  that  shown  as  the  UMAT  curve  of 
figure  12.  Instead  of  abrupt  drops  in  load,  ply  failure  causes  a  reduction  in  the  overall  modulus 
of  the  material.  This  can  be  seen  in  the  slight  deflection  of  the  curve  at  0.35%  and  1%  strain  and 
the  noticeable  deflection  at  3.2%  strain.  To  facilitate  a  one-to-one  validation  of  the  UMAT  and 
the  analytical  code,  simple  modifications  to  the  underlying  FORTRAN  code  of  the  original 
LAM3DNLP  were  made  to  incorporate  this  procedure. 
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Figure  12.  Stress-strain  curves  of  [07907+45°]s  S2-epoxy  laminate  -  LAM3DNLP  output  and 
ABAQUS  UMAT  output. 

3.2.2  Tracking  Progressive  Failure 

Improvements  to  the  tracking  of  failure  progression  were  added  to  the  LAMPATNL  user 
material  subroutine  developed  by  Powers  et  al.  In  LAM3DNLP,  one  material  point  is  simulated, 
but  LAMPATNL  must  account  for  upwards  of  several  thousand  elements  with  multiple 
integration  points.  The  user  material  subroutine  is  called  at  every  integration  point  during  each 
load  increment  within  the  analysis.  During  each  call,  the  UMAT  receives  the  strain  state  of  that 
integration  point  and  determines  the  stiffness  matrix  of  the  material  in  order  for  the  program  to 
calculate  the  stress.  To  accomplish  this,  the  UMAT  must  know  the  failure  history  of  the 
laminate  at  that  integration  point. 

The  failure  history  tracks  which  ply  and  in  which  direction  failure  has  occurred.  This  requires 
the  failure  history  to  be  element,  integration  point,  ply,  and  direction  specific  (figure  13).  With  a 
model  containing  hundreds  of  plies  each  with  six  moduli  (Ei,  E2,  E3,  G12,  G13,  G23)  and 
thousands  of  elements  each  with  multiple  integration  points,  tracking  all  of  this  necessary 
information  can  be  difficult  and  memory  intensive.  It  is  accomplished  by  saving  a  large  array 
that  holds  the  failure  history  and  is  sorted  by  ply  number,  element  number,  and  integration  point 
number. 
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Figure  13.  Tracking  failure  history  is  mode,  ply,  integration  point,  and  element  specific. 

3.2.3  Selection  of  Elements 

Selecting  the  correct  element  is  important  when  simulating  the  behavior  of  complex  structures. 
The  “LAM”  codes  were  developed  to  analyze  thick-section  composite  structures.  To  match  the 
analysis  capability  of  these  programs,  it  is  important  that  LAMPATNL  can  interface  with  3-D 
elements.  LAMPATNL  was  developed  by  Powers  et  al.  to  support  the  use  of  2-D,  solid 
(continuum),  axisymmetrical  four-node  bilinear  (CAX4)  elements.  These  elements  were  first 
supported  because  of  the  need  for  simulating  cylinders  and  other  axisymmetric  structures.  The 
user  material  was  also  able  to  interface  with  elements  with  four-term  stress  and  strain  tensors  of 
the  form  ( <jx , <r  , cr, , crvv )  and  ( ex , sy , c  ,e  ).  To  improve  the  3-D  modeling  capabilities  of 

LAMPATNL,  support  of  eight-node  linear  brick  elements  (C3D8  in  ABAQUS)  was  added.  This 
support  extends  to  any  element  that  has  six-term  stress  and  strain  tensors  in  the  forms 
( cl . <L  > CL > cLy , cr^ , <jyz )  and  ( ex , sy , c_ , , e„ ,£yz). 

It  is  important  to  investigate  which  of  these  supported  elements  are  best  for  accurately  simulating 
thick-section  composite  structures  and  which  provide  easily  interpreted  visualizations  of  key 
outputs.  Modeling  non-axisymmetric  structures  leads  to  the  use  of  the  3-D  eight-node  linear 
brick  elements.  There  are  several  options  available  for  this  node  type,  including  hybrid  with 
both  linear  and  constant  pressure,  the  support  of  incompatible  modes,  and  reduced  integration. 

The  reduced  integration  option  for  this  element  has  several  advantages  and  disadvantages.  Fully 
integrated  C3D8  elements  contain  eight  integration  points  per  element  and  reduced  integration 
C3D8R  elements  contain  only  one  integration  point.  Using  reduced  integration  elements  would 
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result  in  one  eighth  of  the  number  of  calculations  for  the  same  amount  of  elements.  The 
resulting  simulation  would  be  less  accurate,  especially  in  cases  where  reduced  integration 
elements  have  historically  exhibited  problems,  such  as  bending.  Elements  with  reduced 
integration  are  also  prone  to  hourglassing  and  controls  must  be  used  to  prevent  undesirable 
distortions. 


The  drawbacks  of  reduced  integration  elements  are  balanced  by  less  computational  cost  and 
simplified  visualizations.  Figure  14  shows  visualization  of  a  model  using  fully  integrated 
elements  on  the  left  and  reduced  integration  elements  on  the  right.  For  this  case,  the  output  is 
direction  of  the  minimum  Cy  ratio,  MINCIJN,  which  can  be  integer  values  between  zero  and  six. 
Non-integers,  negative  numbers  represented  as  black  areas,  and  number  greater  than  six  (grey 
areas)  do  not  reflect  correct  values.  These  incorrect  values  are  a  product  of  displaying 
extrapolated  values,  derived  from  integration  points,  at  the  nodes  of  full  integration  elements. 

For  this  output,  the  values  can  be  discontinuous  from  one  element  to  the  next,  yet  the  nodal 
extrapolation  assumes  these  values  to  be  continuous.  Using  reduced  integration  elements  avoids 
problems  associated  with  the  visualization  of  extrapolated  nodal  values. 
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Figure  14.  Comparison  of  visualization  of  full  integration  (left)  and  reduced  integration  (right)  elements. 


Careful  consideration  was  taken  to  ensure  all  the  concerns  associated  with  reduced  integration 
elements  are  addressed  and  the  implementation  of  FAMPATNF  in  conjunction  with  these 
elements  was  correct.  Using  UMAT  with  reduced  integration  elements  required  hourglass 
stiffness  controls.  To  ensure  the  hourglass  stiffness  of  the  model  was  set  to  an  appropriate  value, 
the  elastic  strain  energy  was  compared  to  the  hourglass  energy  calculated  by  the  FE  program. 
The  hourglass  energy  should  be  less  than  10%  of  the  elastic  strain  energy  to  ensure  no 
hourglassing  had  occurred.  The  calculation  of  elastic  strain  energy  was  added  to  FAMPATNF 
to  make  this  comparison. 
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3.3  Validation  of  LAMPATNL 

Validating  the  LAMPATNL  UMAT  subroutine  begins  with  modeling  simple  composite  layups 
using  linear  ply  properties  and  no  progressive  failure.  These  single-element  models  are  validated 
with  comparisons  to  the  material  point  response  calculated  by  LAM3D.  Matching  the  stress- 
strain  curve  produced  by  the  UMAT  using  linear  ply  properties  to  that  created  by  LAM3D 
ensures  that  the  user  material  “smears”  the  ply  properties  of  the  laminate  correctly.  While  the 
stress-strain  curve  comparison  is  straightforward,  comparing  the  output  parameters  created  by 
each  program  is  more  complex.  The  parameters  calculated  by  LAM3D  are  the  safety  factor, 
critical  mode,  and  critical  ply.  LAMPATNL  also  produces  these  parameters  and  13  additional 
parameters,  but  there  are  major  differences  in  the  meaning  of  the  parameters  between  the  two. 

For  the  linear  case,  safety  factor  is  defined  as  the  ultimate  load  that  a  laminate  could  sustain, 
based  on  the  allowable  amount  of  progressive  failure,  divided  by  the  applied  loading  on  the 
laminate.  The  safety  factor  is  therefore  a  scalar  multiple  of  the  applied  loading.  The  process  of 
determining  this  value  starts  by  calculating  the  ratio  of  the  strain  allowable  to  ply  strain  for  each 
ply  and  each  direction.  The  lowest  ratio  is  called  the  safety  factor  for  that  iteration  and  the  ply 
and  direction  in  which  the  lowest  ratio  occurs  is  called  the  critical  ply  and  critical  mode, 
respectively.  The  ply  and  direction  are  next  flagged  for  failure  and  the  stiffness  in  the  critical 
direction  of  that  ply  is  discounted.  New  linear  “smeared”  properties  for  the  laminate  are  then 
calculated  and  the  process  repeats.  A  new  safety  factor,  critical  mode,  and  critical  ply  are 
calculated  and  recorded  for  the  next  iteration,  the  plies  are  discounted,  and  “smeared”  properties 
are  updated.  Iterations  continue  until  the  structural  integrity  of  the  laminate  has  been 
compromised.  At  the  end  of  the  analysis,  the  highest  safety  factor  of  all  the  iterations  is 
determined  to  be  last  ply  failure  and  the  mode  and  ply  of  that  highest  safety  factor  are  reported  as 
critical  mode  and  ply  for  the  element. 

The  process  for  determining  the  safety  factor  in  nonlinear,  progressive  failure  analysis  is  not  as 
straightforward.  First,  stress  and  strains  in  nonlinear  analysis  cannot  be  linearly  scaled  to  larger 
values  because  this  analysis  uses  the  incremental  approach  for  determining  stress  from  strain. 
Secondly,  the  progressive  nature  of  the  analysis  requires  that  plies  must  exceed  their  strain 
allowable  to  be  flagged  for  failure.  In  the  linear  analysis,  the  lowest  ratio  of  the  strain  failure 
allowable  to  the  ply  strain  is  designated  for  failure  and  the  ply  is  discounted.  In  nonlinear, 
progressive  failure  analysis,  the  safety  factor  is  defined  as  the  lowest  ratio  of  failure  allowable  to 
ply  strain  at  the  current  level  and  the  critical  ply  and  mode  are  the  ply  and  direction  of  this  lowest 
ratio.  This  determines  which  ply  and  direction  are  closest  to  failure,  not  which  ply  and  direction 
would  cause  last  ply  failure.  Because  of  these  differences,  it  is  not  appropriate  to  directly 
compare  values  of  the  variables  calculated  in  LAMPATNL  to  those  calculated  in  the  LAM3D. 

An  outline  of  the  examples  used  to  validate  the  LAMPATNL  UMAT  subroutine  is  shown  in 
table  3.  Using  linear  ply  properties  and  no  progressive  failure,  two  laminates  are  analyzed  with 
LAMPATNL  and  compared  to  LAM3D  analysis.  With  nonlinear  ply  properties  and  progressive 
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failure,  LAMPATNL  is  validated  against  LAM3DNLP  for  several  materials  and  layups.  For  all 
of  these  examples,  laminate  stress-strain  response  is  compared  between  the  UMAT  and 
analytical  models.  Because  of  the  differences  in  safety  factor  calculations,  this  value  cannot  be 
directly  compared. 

Table  3.  Examples  provided  for  LAMPATNL  validation. 


Lay  Up 

Material 

Load 

(MPa) 

Load  Direction 

Linear/Nonlinear 

Progressive 

Failure 

[0/90/+45] 

S-glass  epoxy 

100 

X  tension 

Linear 

No 

[+30] 

T300/PR-319 

150 

Y  tension 

Linear 

No 

[0] 

S-glass  epoxy 

1700 

1  tension 

Nonlinear 

Yes 

[0] 

S-glass  epoxy 

65 

2  tension 

Nonlinear 

Yes 

[0] 

S-glass  epoxy 

50 

3  tension 

Nonlinear 

Yes 

[0/90/+45] 

S-glass  epoxy 

800 

X  tension 

Nonlinear 

Yes 

[±35] 

IM7/855  1-7 

425 

X  tension 

Nonlinear 

Yes 

[±35] 

IM7/855  1-8 

160 

Y  tension 

Nonlinear 

Yes 

[±35] 

IM7/855  1-9 

85 

Z  tension 

Nonlinear 

Yes 

[±55] 

E-glass/MY750 

125 

X  tension 

Nonlinear 

Yes 

[±55] 

E-glass/MY750 

-200 

X  compression 

Nonlinear 

Yes 

[±55] 

E-glass/MY750 

350 

Y  tension 

Nonlinear 

Yes 

[±55] 

E-glass/MY750 

-225 

Y  compression 

Nonlinear 

Yes 

3.3.1  Validation  of  LAMPATNL  for  Linear  Materials 

In  the  first  example,  an  S-Glass/Epoxy  laminate  with  a  [0°/907±45°]s  layup  is  evaluated  in 
simple  tension  load  of  100  MPa  in  the  X-direction.  The  linear  ply  level  properties  of  this 
material  for  LAM3D  are  shown  in  table  4.  Maximum  strain  is  the  failure  criteria  that  is  used  for 
this  example  and  the  failure  allowables  are  also  shown  in  table  4. 

Table  4.  Linear  ply  properties  of  S-glass/epoxy. 


Ply  Properties 

Strain  Allowables 

E, 

52,000 

£iT 

3.270% 

e2 

19,000 

£iC 

2.210% 

e3 

19,000 

£2T 

0.263% 

G12 

6700 

£2C 

1.500% 

Gi3 

6700 

£3T 

0.263% 

G23 

6700 

£3C 

1.500% 

vi2 

0.30 

y12 

4.000% 

V13 

0.30 

Yi3 

4.000% 

v23 

0.42 

y23 

0.590% 

The  ply  properties  input  into  the  nonlinear  user  material  must  be  linearized,  which  is  achieved  by 
setting  the  stress  asymptote  and  shape  factor  Ramberg-Osgood  parameters  to  large  values. 

Table  5  contains  the  LAMPATNL  input  values  for  the  linearized  material  properties. 
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Additionally,  to  disable  the  progressive  failure  capabilities  of  the  code,  maximum  strain 
allowables  are  set  to  10%. 


Table  5.  Ply  properties  for  linear  S -glass/epoxy. 


S-Glass/Epoxy 

i 

2 

3 

12 

13 

23 

E0  (GPa) 

52 

19 

19 

6.7 

6.7 

6.7 

do  (GPa) 

100 

100 

100 

100 

100 

100 

n 

10 

10 

10 

10 

10 

10 

V 

— 

— 

— 

0.3 

0.3 

0.42 

The  linear  stress-strain  curves  of  both  LAMPATNL  and  the  LAM3D  results  are  shown  in 
figure  15.  The  calculated  stress-strain  responses  from  these  codes  are  coincident  as  expected. 
Because  of  the  reasons  stated  in  the  previous  section,  the  direct  comparison  of  safety  factor, 
critical  mode,  and  critical  ply  are  not  performed. 


Figure  15.  Stress  vs.  strain  curves  in  the  X-direction  for  [0°/90o/+45°]s  S-glass/epoxy  comparing  the 
LAM3D  and  LAMPATNL  codes. 

For  the  next  example,  T300/PR-319  in  a  [±30°]s  layup  is  subjected  to  a  tensile  load  of  150  MPa 
in  the  Y-direction.  The  ply  properties  used  for  the  material  are  listed  in  table  6,  with  the  ao  and  n 
values  not  used  for  input  into  LAM3D.  The  resulting  linear  stress-strain  curves  are  shown  in 
figure  16  and  are  again  coincident. 
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Table  6.  Ply  properties  for  linear  T300/PR-319. 


T  300/PR319 

i 

2 

3 

12 

13 

23 

E0  (GPa) 

129 

5.6 

5.6 

1.33 

1.33 

1.86 

a0  (GPa) 

100 

100 

100 

100 

100 

100 

n 

10 

10 

10 

10 

10 

10 

V 

— 

— 

— 

0.318 

0.318 

0.5 

Figure  16.  Stress  vs.  strain  curves  in  the  Y-direction  for  [+30°]s  T300/PR-319  comparing  the  LAM3D  and 
LAMPATNL  codes. 

3.3.2  Validation  of  LAMPATNL  for  Nonlinear  Materials 

Comparisons  between  the  LAMPATNL  user  material  and  LAM3D  were  made  by  suppressing 
the  nonlinear  and  progressive  failure  features  of  the  user  material.  In  order  to  completely 
validate  its  proper  function,  the  UMAT  model  must  include  these  features.  The  next  step  in 
validating  the  UMAT  is  to  match  results  using  a  single  element  model  with  those  that  are  output 
by  the  analytical  code  LAM3DNLP.  Changes  to  the  post-failure  loading  procedure  for  the 
analytical  code  were  detailed  in  sections  3.2.  These  changes  facilitate  a  one-to-one  comparison 
of  stress-strain  curves,  critical  modes,  and  critical  plies  between  the  user  material  and 
LAM3DNLP. 

The  Ramberg-Osgood  parameters  that  are  used  to  model  the  nonlinear  properties  for  each 
material  used  in  the  examples  are  provided  in  table  7  (30,  35).  The  maximum  strain  failure 
criterion  is  used  for  every  example  and  the  strain  allowables  are  shown  in  table  8  (30,  35). 
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Table  7.  Ramberg-Osgood  parameters. 


Material  and  Its 
Parameters 

Spatial  Directions  for  Constitutive  Modeling 

S-Glass/Epoxy 

1 

2 

3 

12 

13 

23 

E0  (GPa) 

52 

19 

19 

6.7 

6.7 

6.7 

tfo  (GPa) 

100 

0.189 

0.189 

0.071 

0.071 

100 

n 

10 

3.601 

3.601 

2.365 

2.365 

10 

V 

— 

— 

— 

0.3 

0.3 

0.42 

T  300/PR319 

1 

2 

3 

12 

13 

23 

E0  (GPa) 

129 

5.6 

5.6 

1.33 

1.33 

1.86 

CT0  (GPa) 

100 

0.135 

0.135 

0.108 

0.108 

100 

n 

10 

4.728 

4.728 

4.872 

4.872 

10 

V 

— 

— 

— 

0.318 

0.318 

0.5 

IM7/8551-7 

1 

2 

3 

12 

13 

23 

E0  (GPa) 

165 

9.01 

9.01 

5.6 

5.6 

2.8 

ct0  (GPa) 

100 

0.193 

0.193 

0.089 

0.089 

100 

n 

10 

3.604 

3.604 

2.018 

2.018 

10 

V 

— 

— 

— 

0.34 

0.34 

0.5 

E-Glass/MY750 

1 

2 

3 

12 

13 

23 

E0  (GPa) 

45.6 

16.2 

16.2 

6.42 

6.42 

5.79 

<r0  (GPa) 

100 

100 

100 

0.077 

0.077 

100 

n 

10 

10 

10 

1.8 

1.8 

10 

V 

— 

— 

— 

0.278 

0.278 

0.4 

E-Glass/LY556 

1 

2 

3 

12 

13 

23 

E0  (GPa) 

53.5 

17.7 

17.7 

6.36 

6.36 

6.32 

do  (GPa) 

100 

100 

100 

0.076 

0.076 

100 

n 

10 

10 

10 

1.85 

1.85 

10 

V 

— 

— 

— 

0.278 

0.278 

0.4 

AS4/3501-6 

1 

2 

3 

12 

13 

23 

E0  (GPa) 

126 

11 

11 

6.8 

6.8 

3.93 

(T0  (GPa) 

100 

100 

100 

0.097 

0.097 

100 

n 

10 

10 

10 

1.96 

1.96 

10 

V 

— 

— 

— 

0.28 

0.28 

0.4 

Table  8.  Maximum  strain  failure  allowable. 


Material 

Y1T 

(%) 

Y1C 

(%) 

Y2T 

(%) 

Y2C 

(%) 

Y3T 

(%) 

Y3C 

(%) 

Y23 

(%) 

Y13 

(%) 

Y12 

(%) 

S-glas  s/epoxy 

3.27 

-2.21 

0.33 

-1.5 

0.263 

-1.5 

0.59 

4 

4 

T  300/PR319 

1.07 

-0.74 

0.43 

-2.8 

0.43 

-2.8 

1.5 

8.6 

8.6 

IM7/8551-7 

1.551 

-0.96 

0.87 

-3.2 

0.755 

-3.2 

2.1 

5 

5 

E-glass/MY750 

2.81 

-1.75 

0.25 

-1.2 

0.25 

-1.2 

4 

4 

4 

E-glass/LY556 

2.13 

-1.07 

0.2 

-0.64 

0.2 

-0.64 

3.8 

3.8 

3.8 

AS4/3501-6 

1.38 

-1.18 

0.44 

-2 

0.44 

-2 

2 

2 

2 

Using  a  single-element  model  for  each  example,  the  laminate  is  loaded  in  the  X-direction  until 
failure  has  occurred.  This  process  is  similarly  repeated  for  the  Y-  and  Z-directions.  The  first 
comparison  between  the  LAMPATNL  user  material  and  LAM3DNLP  is  a  simple  unidirectional 
layup  of  S- Glass/Epoxy.  As  with  all  of  the  materials  tested,  this  material  is  linear  in  the  fiber 
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direction.  The  tensile  failure  strains  in  the  2-  and  3-directions  for  this  material  are  relatively 
small,  which  results  in  the  nonlinear  responses  of  the  material  in  those  directions  appearing  to  be 
linear.  Figure  17  shows  the  stress-strain  response  of  the  unidirectional  material  in  the  1- 
direction.  As  expected,  the  UMAT  curve  is  coincident  with  the  analytical  model  and  both  reach 
failure  at  3.27%  strain. 


Figure  17.  Stress  vs.  strain  curves  in  the  1-direction  for  [0°]  S-glass/epoxy  comparing  the  UMAT  and 
LAM3DNLP  codes. 

The  stress-strain  response  in  the  2-direction  is  shown  in  figure  18.  Failure  in  this  direction 
occurs  at  0.33%  strain  and  this  is  reflected  in  the  results.  Figure  19  shows  the  response  of  this 
material  in  the  3-direction,  with  failure  occurring  at  0.26%  strain  in  both  the  LAMPATNL  and 
LAM3DNLP. 


34 


Strain  (%) 


Figure  18.  Stress  vs.  strain  curves  in  the  2-direction  for  [0°]  S-glass/epoxy  comparing  the  UMAT  and 
LAM3DNLP  codes. 


Strain  (%) 


Figure  19.  Stress  vs.  strain  curves  in  the  3-direction  for  [0°]  S-glass/epoxy  comparing  the  UMAT  and 
LAM3DNLP  codes. 
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The  next  example  is  an  S-Glass/Epoxy  laminate  in  a  [07907+45 °]s,  or  quasi-isotropic,  layup.  In 
this  example,  the  response  to  load  applied  in  the  X-direction  is  the  same  as  response  to  load 
applied  in  the  Y-direction,  so  only  one  plot  is  shown  for  these  two  load  cases.  The  nonlinear 
material  response  and  progressive  failure  of  the  laminate  is  shown  in  figure  20.  Transverse 
tensile  failure  in  the  90°  plies  occurs  in  the  laminate  at  0.35%  strain  and  in  the  ±45°  plies  at  1% 
strain.  These  failures  produce  slight  deflections  in  the  stress-strain  response.  Some  additional 
failure  has  little  effect  on  the  overall  stiffness  of  the  laminate  in  this  direction.  Longitudinal 
tensile  failure  in  the  0°  plies  occurred  at  3.3%  strain  that  caused  a  dramatic  decrease  in  the 
stiffness  of  the  laminate.  The  final  failure  recorded  for  this  case  is  longitudinal  compressive 
failure  in  the  90°  plies.  Both  the  analytical  code  and  the  user  material  yield  the  same  results  for 
this  test  case,  continuing  to  validate  the  correct  function  of  the  user  material. 


Strain  (%) 


Figure  20.  Stress  vs.  strain  curves  in  the  X-direction  for  [0°/90°/+45°]s  S-glass/epoxy  comparing  the  UMAT 
and  LAM3DNLP  codes. 

The  out-of-plane  stiffness  of  the  laminate  is  not  significantly  affected  by  its  layup  and  there  is  no 
progressive  failure  in  this  direction  because  first  ply  failure  is  catastrophic.  The  stress-strain 
response  in  the  Z-direction  is  the  same  as  that  shown  in  figure  19  because  the  material  used  for 
this  laminate  is  the  same  used  in  the  first  example. 

A  laminate  comprised  of  IM7/8551-7  in  a  [±35°]s  layup  is  used  in  the  next  validation  test  case. 

In  this  composite  layup,  the  nonlinear  transverse  and  shearing  directions  are  a  contributing  factor 
to  the  response  in  X-direction  loading.  The  stress-strain  curves  in  figure  21  do  not  show  abrupt 
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discontinuities  or  changes  in  stiffness  outside  those  normally  associated  with  nonlinear  materials, 
suggesting  that  the  failures  did  not  significantly  affect  the  response  of  the  laminate  in  this 
direction. 


Strain  (%) 


Figure  21.  Stress  vs.  strain  curves  in  the  X-direction  for  [+35°]s  IM7/8551-7  comparing  the  UMAT  and 
LAM3DNLP  codes. 

Subjecting  the  same  laminate  to  a  load  in  the  Y-direction  yields  a  response  that  is  both  highly 
nonlinear  and  significantly  affected  by  progressive  ply  failure.  In  figure  22,  the  stress-strain 
response  deviates  from  the  nonlinear  behavior  to  a  more  compliant  response  at  approximately 
1.7%  strain.  This  failure  corresponds  to  transverse  tensile  failure  in  the  ±35°  plies.  Ultimate 
failure  in  this  laminate  is  in-plane  shearing  in  all  the  plies.  The  results  from  the  user  material 
once  again  mirror  the  results  obtained  by  LAM3DNLP. 
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Figure  22.  Stress  vs.  strain  curves  in  the  Y-direction  for  [+35°]s  IM7/8551-7  comparing  the  UMAT  and 
LAM3DNLP  codes. 

The  response  of  this  laminate  in  the  Z-direction  (figure  23)  is  linear  until  failure.  Once  again,  the 
small  tensile  failure  strain  of  the  material  limits  its  response  to  the  linear  regime  of  the  nonlinear 
curve.  Both  the  user  material  and  LAM3DNLP  responses  continue  to  match. 

The  final  example  for  validation  is  an  E-Glass/MY750  composite  in  a  [±55°]s  layup.  This 
example  tests  the  X-  and  Y-direction  response  of  the  laminate  in  both  tension  and  compression. 
The  material  properties,  layup,  and  load  scenario  are  the  same  as  those  studied  in  the  first  World 
Wide  Failure  Exercise  (WWFE-I  (see  test  case  9)  (30).  The  details  of  this  exercise  are  discussed 
in  the  Background  of  the  “LAM”  codes  section  of  section  2. 


38 


Figure  23.  Stress  vs.  strain  curves  in  the  Z-direction  for  [+35°]s  IM7/8551-7  comparing  the  UMAT  and 
LAM3DNLP  codes. 

The  first  validation  for  this  laminate  is  tension  in  the  X-direction.  In  figure  24,  the  laminate  has 
an  almost  purely  linear  response  up  to  around  0.5%  strain  then  a  noticeable  shift  to  nonlinear 
response  after.  This  shift  is  caused  by  transverse  tensile  failure  in  all  the  plies  of  the  laminate. 
The  coincidence  of  the  curves  further  validates  the  results  of  the  UMAT. 
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Figure  24.  Stress  vs.  strain  curves  in  X-direction  tension  for  [+55°]s  E-Glass/MY750  comparing  the  UMAT 
and  LAM3DNLP  codes. 

The  second  validation  is  compression  in  the  X-direction,  shown  in  figure  25.  Like  previous 
examples,  the  LAM3DNLP  results  and  the  UMAT  results  match.  Transverse  compressive 
failure  in  the  plies  around  2.5%  strain  causes  the  laminate  to  become  compliant.  Ultimate  failure 
in  this  test  case  is  a  result  of  shear  failure  in  the  plies. 
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Figure  25.  Stress  vs.  strain  curves  in  X-direction  compression  for  [+55°]s  E-Glass/MY750  comparing  the 
UMAT  and  LAM3DNLP  codes. 

Tension  in  the  Y-direction  is  next  investigated  for  validation  and  shown  in  figure  26.  The 
response  of  the  laminate  is  nonlinear  with  no  noticeable  stiffness  changes  that  can  be  attributed 
to  progressive  failure.  In-plane  shear  failure  in  this  laminate  around  2.25%  strain  does  not 
significantly  affect  the  response  in  this  direction.  Ultimate  failure  in  this  test  case  is  a  result  of 
transverse  compressive  failure. 
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Figure  26.  Stress  vs.  strain  curves  in  Y-direction  tension  for  [+55°] s  E-Glass/MY750  comparing  the 
UMAT  and  LAM3DNLP  codes . 

The  final  example  in  figure  27  is  compression  in  the  Y-direction.  A  change  in  stiffness  around 
1.25%  strain  is  a  result  of  transverse  tensile  failure  in  the  plies.  This  failure  has  a  small  effect  on 
the  response  of  the  laminate  in  this  direction.  In-plane  shear  failure  in  all  plies  is  the  ultimate 
failure  in  this  case.  Similar  to  all  other  validation  examples,  the  UMAT  stress-strain  curves 
match  those  produced  by  the  analytical  LAM3DNLP  code. 
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Figure  27.  Stress  vs.  strain  curves  in  Y-direction  compression  for  [+55°]s  E-Glass/MY750  comparing 
the  UMAT  and  LAM3DNLP  codes. 


3.4  Validation  Summary 

The  comparison  of  the  LAMPATNL  user  material  output  using  a  single-element  model  to 
LAM3D  for  linear  properties  provided  the  first  step  towards  validating  the  proper  function  of  the 
subroutine.  Four  example  cases  were  shown  to  validate  the  nonlinear  and  progressive  failure 
aspects  of  the  user  material  function  correctly  by  generating  the  same  output  as  LAM3DNLP. 
The  validated  LAMPATNL  UMAT  can  now  be  used  in  the  design  of  composite  structures.  The 
procedure  of  incorporating  the  UMAT  into  the  design  process  of  thick  section  composites  is 
documented  in  section  4. 


4.  Case  Studies 


LAM3DNLP  and  LAMPATNL  were  developed  to  analyze  symmetric,  thick-section  composite 
laminates  and  structures  with  ply  level  material  nonlinearity  and  failure.  In  this  section,  two  case 
studies  are  presented  and  examined  with  the  LAMPATNL  UMAT.  The  first  case  study  is  a 
simple  block  in  compression  that  is  partially  fixed  on  one  edge,  inducing  in-plane  shear  stresses. 
The  second  case  study  is  an  open  hole  sample  subjected  to  simultaneous  loading  in  the  X-  and  Y- 
directions.  The  open  hole  sample  is  commonly  used  to  check  the  effectiveness  of  composite 
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failure  models  (14,  15,  17-21,  23,  25,  and  26).  An  example  illustrating  the  improvements  and 
contributions  to  the  LAMPATNL  design  method  over  its  original  version  is  also  provided. 

4.1  LAMPATNL  Output  Parameters 

A  total  of  16  parameters  are  output  by  the  UMAT,  so  filtering  and  processing  this  data  are 
important  parts  of  the  design  process.  Several  useful  output  parameters  are  calculated  for  each 
element  in  the  model  in  addition  to  the  safety  factor,  critical  mode,  and  critical  ply  parameters 
produced  by  the  original  LAMPATNL  user  material  (36).  A  full  list  and  brief  description  of 
each  parameter  is  provided  in  table  9.  These  history  dependent  parameters  are  calculated  and 
stored  throughout  the  incremental  solution. 

Table  9.  List  of  state  dependent  parameters. 


Abaqus  Designation 

Parameter  Name 

Description 

SDV1 

CXXR 

Cxx  ratio,  X  stiffness 

SDV2 

CYYR 

Cxx  ratio,  Y  stiffness 

SDV3 

CZZR 

Cxx  ratio,  Z  stiffness 

SDV4 

CYZR 

Cxx  ratio,  YZ  stiffness 

SDV5 

CXZR 

Cxx  ratio,  XZ  stiffness 

SDV6 

CXYR 

Cxx  ratio,  XY  stiffness 

SDV7 

NOF 

Number  of  failures 

SDV8 

CMODE 

Current  critical  failure  mode 

SDV9 

FMODE 

Previous  failed  mode 

SDV10 

CPLY 

Current  critical  ply  number 

SDV11 

FPLY 

Previous  failed  ply  number 

SDV12 

SF 

Safety  factor 

SDV13 

SUMCIJ 

Sum  of  SDV1-SDV6 

SDV14 

PRODCIJ 

Product  of  SDV1-SDV6 

SDV15 

MINCIJ 

Minimum  value  of  SDV1-SDV6 

SDV16 

MINCIJN 

Cn  number  (1-6)  ofSDV15 

The  first  six  parameters  are  stiffness  ratios  evaluated  directly  from  the  elemental  (global) 
stiffness  matrix.  At  the  beginning  of  the  FE  analysis,  the  values  on  the  diagonal  of  the  stiffness 
matrix,  Qj,  are  stored.  These  values  change  throughout  the  analysis  due  to  material  nonlinearity 
and  progressive  ply  failure.  The  stiffness  matrix  is  updated  and  the  current  values  occupying  the 
diagonal  are  divided  by  the  original  values  to  determine  the  stiffness  ratio  for  that  direction. 
Equation  38  shows  the  calculation  of  the  Cxx  stiffness  ratio  CXXR,  with  the  calculation  of  the 
five  other  ratios  following  similarly. 

CXXR  =  %-  (38) 

Cxx 
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In  equation  38,  the  c  term  is  the  current  value  and  the  i  term  is  the  initial  value.  For  this  analysis, 
Cxx  refers  to  value  in  the  (1,  1)  position  of  the  global  Qj  stiffness  matrix,  Cyy  in  the  (2,  2),  Czz 
in  the  (3,  3),  Cyz  in  the  (4,  4),  Cxz  in  the  (5,  5),  and  Cxy  in  the  (6,  6).  The  values  of  the  stiffness 
ratios  range  from  1  to  0,  where  a  value  of  1  indicates  no  stiffness  loss  from  the  original  stiffness 
and  0  indicates  complete  loss  of  stiffness  of  the  element  in  the  direction. 

The  parameter  for  the  number  of  failures,  NOF  or  SDV7,  records  the  number  of  ply  failures  that 
have  occurred  in  the  element.  For  example,  if  a  [0790°] s  laminate  had  transverse  tensile  failure 
in  the  90°  plies  and  longitudinal  tensile  failure  in  the  0°  plies  the  number  of  failures  would  be  4. 
This  value  is  a  result  of  symmetry  with  a  single  failure  recorded  in  the  two  90°  plies  and  another 
single  failure  recorded  in  the  two  0°  plies.  If  failure  has  not  occurred  in  the  element,  NOF  will 
be  at  zero.  This  parameter  does  not  indicate  the  effect  of  ply  failure  on  the  stiffness  of  the 
sample  but  is  a  count  of  the  number  of  failure  allowables  exceeded  in  the  element  up  to  this 
point. 


The  critical  mode  CMODE,  previous  failed  mode  FMODE,  critical  ply  CPLY,  previous  failed 
ply  FPLY,  and  safety  factor  SF  are  important  parameters  used  to  determine  the  progressive  ply 
failure  history  in  LAMPATNL.  The  safety  factor  SF  is  defined  as  the  ratio  of  the  principal 
material  direction  strain  allowable  to  the  current  ply  level  strain,  as  shown  below: 


YIT 

sf ;  = 

£l  if  8i>0 

(39) 

Y1C 

SF2  =  — 

£l  if  8i<0 

(40) 

Y2T 

SF3  = 

£l  if  £2>0 

(41) 

Y2C 

sf4  =  YZL 

£2  if  82<0 

(42) 

Y3T 

sf5=  — 

£i  if  £3>0 

(43) 

Y3C 

sf6  =  — 

if  £3<0 

(44) 

(45) 
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(46) 


SFs 


Y 13 

^5 


sf9  = 


(47) 


where  Y1T,  Y1C,  Y2T,  Y2C,  Y3T,  Y3C,  Y23,  Y13,  and  Y12  are  the  maximum  strain 
allowables.  These  calculations  are  done  for  all  N  plies  in  the  laminate.  The  lowest  value  out  of 
the  9N  safety  factors  is  recorded  and  displayed  as  the  SF  output  parameter.  The  ply  with  the 
lowest  SF  becomes  the  critical  ply  CPLY  and  the  mode  in  which  this  SF  has  occurred  is  recorded 
as  the  critical  mode  CMODE.  The  values  of  CPLY  range  from  1  to  N.  The  values  of  CMODE 
correspond  to  the  number  of  the  lowest  SF,  as  defined  in  equations  39-47,  and  range  from  1  to  9. 
Critical  or  failure  mode  1  refers  to  1 -direction  fiber  tension,  mode  2  is  1 -direction  fiber 
compression,  mode  3  is  2-direction  matrix  tension,  mode  4  is  2-direction  matrix  compression, 
mode  5  is  3-direction  matrix  tension,  mode  6  is  3-direction  matrix  compression,  mode  7  is 
interlaminar  (23)  shear,  mode  8  is  interlaminar  (13)  shear,  and  mode  9  is  in-plane  (12)  shear. 

The  failure  modes  refer  to  the  ply-level  coordinate  system  that  is  related  to  the  laminate  and 
global  systems  as  shown  in  figure  28. 


Figure  28.  Relation  of  global  coordinates  to  laminate  and  ply  coordinates. 

When  a  strain  in  a  ply  has  reached  its  allowable  (SF=1),  the  failed  ply  and  direction  are 
discounted  and  excluded  from  the  SF  calculation.  The  critical  ply  at  this  point  is  recorded  as  the 
previous  failed  ply  number  FPLY  and  the  critical  mode  becomes  the  previous  failed  mode 
FMODE.  These  values  are  initially  at  zero  and  will  take  on  the  range  of  values  of  CPLY  and 
CMODE  if  ply  failure  occurs  in  the  element.  If  failure  is  recorded,  the  FPLY  and  FMODE 
parameters  will  portray  the  progression  of  failure  in  the  analysis. 
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The  SUMCIJ  parameter  is  simply  the  sum  of  the  CXXR,  CYYR,  CZZR,  CYZR,  CXZR,  and 
CXYR  stiffness  ratio  parameters  and  can  range  from  6  to  0.  The  PRODCIJ  parameter  is  the 
product  of  these  stiffness  ratios  and  ranges  from  1  to  0.  The  MINCU  parameter  records  the 
minimum  value  of  the  six  stiffness  ratios  and  also  ranges  from  1  to  0.  The  MINCIJN  parameter 
indicates  which  stiffness  ratio  is  recorded  as  the  minimum  for  the  element.  If  CXXR  is  the 
minimum  stiffness  ratio,  MINCIJN  is  1,  for  CYYR  it  is  2,  for  CZZR  it  is  3,  for  CYZR  it  is  4,  for 
CXZR  it  is  5,  and  for  CXYR  it  is  6.  If  the  minimum  stiffness  ratio  for  the  element  is  greater  than 
80%,  no  value  for  MINCUN  is  recorded  and  this  parameter  is  set  to  zero.  Some  output 
parameters  are  not  used  in  the  following  case  studies  but  may  have  more  important  uses  in 
different  analyses. 

4.2  Design  Methodology 

A  standardized  process  for  designing  composite  structures  with  the  LAMPATNL  user  material 
has  been  developed.  The  output  parameters  produced  by  the  user  material  are  subsequently  used 
to  analyze  the  nonlinear  response  and  progressive  failure  of  the  composite  structure.  The 
methodology  for  designing  a  structure  using  LAMPATNL  is  outlined  in  figure  29. 


Figure  29.  Design  methodology  flowchart. 
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The  largest  reduction  in  original  stiffness,  identified  by  the  lowest  Qj  ratio  MINCU  at  the  end  of 
analysis  and  its  corresponding  direction  MINCUN,  are  examined  first.  The  next  step  is  to 
determine  from  the  NOF  parameter,  which  records  the  number  of  failures,  whether  the  reduction 
in  stiffness  has  been  caused  by  material  nonlinearity  or  ply  failures.  If  there  are  ply  failures  in 
the  areas  of  concern,  the  complete  failure  progression  of  the  structure  is  constructed  using  the 
failed  mode  FMODE  and  failed  ply  FPLY.  If  the  areas  of  low  stiffness  ratio  are  caused  by 
nonlinear  effects,  the  critical  ply  CPLY  and  mode  CMODE  at  the  final  increment  become  the 
basis  for  changes  to  the  layup.  The  designer  would  use  this  information,  critical  ply  and  mode  or 
failure  progression,  to  make  changes  to  the  design  of  the  structure  to  achieve  the  objectives  of 
the  analysis.  Multiple  iterations  of  this  design  methodology  may  be  needed  to  achieve  these 
objectives.  The  following  case  studies  help  to  demonstrate  this  design  methodology. 

4.3  Original  LAMPATNL  Design  Methodology 

An  example  comparing  the  original  outputs  of  LAMPATNL  (safety  factor,  critical  mode,  and 
critical  ply)  to  the  new  stiffness  ratios  outputs  is  provided  in  this  section.  This  example 
illustrates  new  contributions  that  the  stiffness  ratio  parameters  bring  to  visualizing  the  critical 
design  information  of  the  structure. 

In  this  example,  a  laminate  of  a  [0°/90°]s  layup  of  IM7/8551-7  is  tested  as  a  beam  in  bending. 
The  Ramberg-Osgood  parameters  and  maximum  strain  allowables  for  this  material  are  shown  in 
tables  6  and  7,  respectively.  The  sample  is  constrained  in  X-  and  Z-directions  along  the  bottom 
right  edge  and  the  Z-direction  along  the  bottom  left  edge  as  shown  in  figure  30.  Also  in  figure 
30,  a  displacement  in  the  Z-direction  along  a  portion  of  the  top  face  is  applied  to  introduce  a 
bending  load  into  the  sample. 


Figure  30.  Diagram  of  a  bending  sample. 
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The  outputs  that  were  originally  available  in  LAMPATNL  analysis  were  the  safety  factor  SF, 
critical  mode  CMODE,  and  critical  ply  CPLY.  When  a  strain  in  a  ply  has  reached  its  allowable, 
the  ply  fails.  This  ply  and  direction  are  discounted  and  excluded  from  the  SF  calculation.  The 
recorded  value  of  the  SF  therefore  decreases  to  one  before  jumping  to  a  different  value,  as  shown 
for  the  loading  of  an  arbitrary  element  in  figure  31. 


Displacement  (cm) 


Figure  31.  Plot  of  safety  factor  vs.  displacement. 

At  the  end  of  the  simulation  for  this  example,  the  element  shown  in  figure  31  displays  a  SF  of 
4.53  yet  there  have  been  at  least  three  ply  failures  in  this  element.  Information  on  the  failure 
history  of  the  laminate  in  this  element  cannot  be  determined  from  this  final  value  of  SF. 

Figure  32  shows  the  contour  plots  of  SF  from  the  beginning  to  the  end  of  the  analysis.  In  these 
plots,  blue  areas  have  a  SF  close  to  1,  red  areas  indicate  a  SF  of  5,  and  grey  areas  are  anything 
greater  than  5.  From  the  plot  of  the  final  increment  (last  plot  in  the  series),  it  can  be  determined 
that  the  blue  areas  are  close  to  failure,  but  there  is  no  information  on  which  areas  and  plies  have 
failed  nor  is  there  information  about  how  these  failures  have  affected  the  stiffness  of  the 
structure.  The  complex  and  abrupt  fluctuations  in  the  SF  and  the  inability  to  determine  the 
failure  history  prevent  the  original  outputs  of  LAMPATNL  from  effectively  providing  the 
information  necessary  to  make  changes  to  the  composite  design. 
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Figure  32.  Plot  of  progression  of  SF. 

In  this  work,  new  output  parameters  from  LAMPATNL  are  proposed  that  include  the  structural 
stiffness  ratios  of  the  laminate.  These  parameters  are  detailed  in  the  LAMPATNL  output 
parameters  section  of  this  chapter.  Figure  33  shows  the  results  of  the  MINCIJ  parameter  in  the 
same  arbitrary  element  documented  in  figure  31.  This  plot  shows  that  there  was  catastrophic 
stiffness  loss  in  at  least  one  direction  for  this  element.  At  the  end  of  the  analysis,  the  value  of 
this  parameter  is  0.007,  indicating  that  complete  stiffness  loss  has  occurred  in  at  least  one 
direction  in  this  element. 
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Figure  33.  Plot  of  minimum  stiffness  ratio  (MINCIJ)  vs.  displacement. 

A  progression  plot  of  the  MINCIJ  parameter  is  shown  in  figure  34.  In  this  plot,  areas  that  are 
blue  correspond  to  a  zero  stiffness  ratio  (complete  stiffness  loss)  and  areas  that  are  red 
correspond  to  no  stiffness  loss.  By  examining  the  contour  plot  of  the  last  increment,  it  is  easy  to 
determine  which  areas  have  severe  stiffness  loss  and  which  areas  do  not.  Because  of  the  simple 
nature  of  the  stiffness  ratio  measure,  it  is  easy  to  determine  from  the  progression  which  areas 
have  failed  and  in  which  order  they  have  failed.  This  contrasts  from  the  SF  plots  shown  in 
figure  32,  which  do  not  efficiently  visualize  the  failure  history  of  the  part. 
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Figure  34.  Plot  of  progression  of  minimum  stiffness  ratio  (MINCIJ). 

4.4  Case  Study  No.  1 
4.4.1  Description 

The  first  case  study  that  is  examined  with  the  LAMPATNL  is  a  compressive  shear  sample.  This 
sample  was  created  to  represent  a  composite  material  that  is  bonded  to  a  much  stiffer  material. 
The  sample,  shown  in  figure  35,  is  rectangular  with  an  aspect  ratio  of  2.5  that  is  fixed  on  one  half 
of  one  side.  A  compressive  load  in  the  X-direction  is  applied  to  one  end  of  the  sample  and  a 
boundary  condition  on  the  opposing  end  prevents  the  sample  from  moving  in  the  X-direction. 
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Figure  35.  Diagram  of  shear  sample. 

The  first  design  tested  is  a  [0790°]s  layup  of  AS4/3501-6.  The  Ramberg-Osgood  parameters  and 
maximum  strain  allowables  for  this  material  are  shown  in  tables  6  and  7,  respectively.  A 
uniform  pressure  of  800  MPa  is  applied  to  the  top  end  of  the  sample,  causing  compressive 
stresses  and  a  shear  stress  concentration  close  to  the  right  fixed  boundary  condition. 

The  sample  is  modeled  using  a  mesh  of  6448  3-D  eight-node  linear  brick  elements  with  reduced 
integration,  type  C3D8R  in  ABAQUS.  The  mesh  contains  8505  nodes  each  with  three  degrees 
of  freedom,  resulting  in  25515  total  degrees  of  freedom.  For  each  case  study,  the  mesh  was 
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varied  in  both  density  and  structure  to  ensure  the  results  of  the  simulation  did  not  change 
dramatically. 

The  objective  of  this  case  study  is  to  limit  the  displacement  of  point  A  (see  figure  35)  in  the  X- 
direction  to  0.05  m.  Limiting  the  shear  and  longitudinal  failure  is  important  to  the  performance 
of  this  sample,  so  reducing  areas  of  stiffness  loss  in  these  directions  is  also  an  objective.  Based 
upon  the  performance  of  the  first  layup,  modifications  to  the  laminate  may  be  needed  to  achieve 
the  objective.  Any  change  to  the  laminate  must  maintain  the  same  material,  weight,  and 
thickness  from  the  original  layup.  These  changes  are  made  based  on  information  gained  from 
the  process  outlined  in  the  design  methodology. 

4.4.2  Results  of  Original  Design 

Following  the  procedure  outlined  in  the  design  methodology,  the  first  step  is  to  investigate  the 
lowest  ratio  of  current  stiffness  to  original  stiffness  at  the  end  of  the  analysis.  The  contour  plot 
of  the  LAMPATNL  output  parameter  MINCIJ  is  shown  in  figure  36.  In  this  plot,  areas  of 
concern  are  those  that  are  blue,  corresponding  to  a  zero  stiffness  ratio,  and  areas  of  least  concern 
are  red,  corresponding  to  no  stiffness  loss.  This  color  scale  holds  for  all  plots  of  the  MINCU, 
CXXR,  CYYR,  CZZR,  and  CXYR  shown  in  this  case  study. 

A  significant  area  of  the  structure  has  at  least  one  of  its  stiffness  ratios  near  zero.  The  next  step 
in  this  methodology  is  to  determine  which  direction  the  stiffness  ratios  are  at  or  near  zero.  A  plot 
of  the  Cij  number  (1-6),  MINCIJN,  indicates  the  critical  direction  of  the  lowest  stiffness  ratio  and 
is  shown  in  figure  37. 
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Figure  36.  Contour  plot  of  minimum  C,j  stiffness  ratio,  MINCIJ. 
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Figure  37.  Contour  plot  of  the  direction  of  minimum  Qj  ratio,  MINCIJN. 


The  area  in  the  bottom  right  portion  of  the  structure  at  the  fixed  half-face  has  Cxy  (red  area  of 
MINCIJN=6)  as  the  lowest  stiffness  ratio,  corresponding  to  the  XY  shear  stiffness.  The  areas  in 
green  also  have  a  stiffness  ratio  close  to  zero  and  correspond  to  a  value  of  3  in  this  plot.  This 
result  shows  that  there  is  a  loss  of  stiffness  in  the  Czz  stiffness  ratio,  or  global  Z-direction.  This 
is  not  a  result  of  loading  in  the  Z-direction,  but  a  result  of  ply  strains  caused  by  Poisson’s  effects 
exceeding  the  low  through-thickness  tensile  failure  allowable. 

In  order  to  further  investigate  the  critical  losses  in  stiffness,  a  closer  examination  of  each 
individual  stiffness  ratio  is  required.  Upon  inspection,  the  CXXR,  CYYR,  CZZR,  and  CXYR 
stiffness  ratio  parameters  are  low  in  the  areas  of  interest.  The  contour  plot  of  the  CZZR  ratio  in 
figure  38  confirms  what  is  shown  in  figure  37,  that  a  large  part  of  the  sample  has  lost  stiffness  in 
the  Z-direction.  The  plot  of  the  Cxy  stiffness  ratio,  shown  in  figure  39,  also  confirms  what  is 
shown  in  figure  37,  that  the  XY  shear  stiffness  is  reduced  near  the  fixed  face  of  the  sample. 
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Figure  38.  Plot  of  Czz  stiffness  ratio,  CZZR,  case  study  no.  1. 
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Figure  39.  Plot  of  CXy  stiffness  ratio,  CXYR,  case  study  no.  1. 


A  closer  examination  of  the  plot  of  the  Cyy  stiffness  ratio  shows  that  the  upper  portion  of  the 
fixed  face  also  experiences  a  loss  of  stiffness  in  the  Y-direction  as  shown  in  figure  40.  The 
limitation  of  the  contour  plot  in  figure  37  is  that  only  one  value  for  each  element  can  be 
displayed.  In  the  upper  portion  of  the  fixed  face,  the  loss  of  stiffness  occurs  in  both  the  Y-  and 
XY-directions  although  only  the  shear  direction  is  indicated  in  this  area  in  figure  37. 
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Figure  40.  Plot  of  Cyy  stiffness  ratio,  CYYR,  case  study  no.  1. 


A  similar  scenario  occurs  in  the  areas  of  figure  37  that  indicate  Z-direction  failure.  Figure  41 
shows  the  plot  of  the  Cxx  stiffness  ratio,  in  which  a  large  part  of  the  area  that  had  stiffness  loss  in 
the  Z-direction  also  has  a  loss  of  X-direction  stiffness.  Because  the  sample  is  loaded  in  this 
direction,  this  failure  is  important  to  the  response  of  the  structure. 
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Figure  41.  Plot  of  CXx  stiffness  ratio,  CXXR,  case  study  no.  1. 


The  next  step  of  the  design  methodology  is  to  determine  if  there  are  any  ply  failures  recorded. 
Because  a  large  part  of  the  sample  has  a  stiffness  ratio  at  or  near  zero,  ply  failures  are  to  be 
expected.  A  plot  of  the  number  of  failures  that  have  occurred  is  shown  in  figure  42.  The  upper 
portion  of  the  fixed  face  has  the  most  ply  failures.  This  plot  also  shows  areas  that  have  a  loss  of 
stiffness  in  the  X-direction  are  differentiated  from  those  with  only  Z-direction  failure,  shown  as 
the  areas  of  orange  and  green,  respectively. 
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Figure  42.  Contour  plot  of  NOF,  case  study  no.  1. 


Continuing  in  the  design  methodology,  the  ply-level  failure  history  of  the  sample  is  investigated. 
Ply-level  failure  behavior  is  the  cause  of  the  global  stiffness  loss  shown  in  the  stiffness  ratio  plots 
in  figures  38-41.  Figure  43  displays  the  progression  of  the  modes  that  are  failing  throughout  the 
response  of  the  structure.  The  failure  modes  shown  in  this  figure  refer  to  the  ply-level  coordinate 
system  as  described  in  section  4.1. 

In  figure  43,  only  the  mode  of  failure  for  the  last  failed  ply  is  displayed.  Areas  in  dark  blue 
(FMODE=0)  have  not  failed,  red  (FMODE=9)  indicates  a  failure  mode  of  12  shear,  green  areas 
(FMODE=5)  have  tensile  failure  in  the  3-direction,  the  areas  in  pale  green  (FMODE=3)  are  a  2- 
direction  tensile  failure,  and  the  light  blue  areas  (FMODE=2)  indicate  a  1 -direction  compressive 
failure.  It  is  important  to  understand  which  of  the  four  plies  correspond  to  the  failure  modes  of 
figure  43.  A  progression  plot  of  the  last  failed  ply  is  shown  in  figure  44. 


61 


ON 


SDV9 


i 


I 


+9.000e+00 
+8.625e+00 
+8.  250e+00 
+7.875e+00 
+7.500e+00 
+7.125e+00 
+6.750e+00 
+6.375e+00 
+6. 000e+00 
+5.625e+00 
+5.250e+00 
+4.875e+00 
+4.500e+00 
+4. 125e+00 
+3.750e+00 
+3.375e+00 
+3. 000e+00 
+2.625e+00 
+2.250e+00 
+1.875e+00 
+1.500e+00 
+1.125e+00 
+7.500e-01 
+3.750e-01 
+0. 000e+00 


Figure  43.  Plot  of  progression  of  last  failed  mode,  FMODE,  case  study  no.  1. 
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Figure  44.  Plot  of  progression  of  last  failed  ply,  FPLY,  case  study  no.  1. 


In  figure  44,  the  blue  areas  indicate  that  no  failure  has  occurred,  red  areas  (FPLY=4)  indicate  that 
the  failure  shown  in  figure  43  occurred  in  the  0°  plies,  and  the  light  orange  area  (FPLY=3) 
indicates  that  the  failure  occurred  in  the  90°  plies.  In  order  to  describe  the  progression  of  failure 
in  the  sample,  four  regions  are  identified  and  highlighted  on  the  sample  in  figure  45.  The 
regions,  shown  in  figure  45  on  a  plot  of  FMODE  at  the  end  of  the  analysis,  are  used  to  refer  to 
areas  in  the  progression  plots  of  figure  43a-e  and  figure  44a-e. 


Figure  45.  Diagram  of  regions,  case  study  no.  1. 

There  are  three  phases  of  failure  that  occur  in  this  sample.  The  first  phase  is  in-plane  shear 
failure  in  all  of  the  plies  that  propagates  from  region  1,  the  top  edge  of  the  fixed  face,  downwards 
into  region  2.  This  phase  is  displayed  in  plots  a  and  b  of  figures  43  and  44.  The  next  phase 
includes  tensile  failure  in  the  3-direction  in  all  plies  that  starts  at  the  top  of  the  fixed  face 
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(region  1)  and  continues  up  the  right  edge  of  the  sample  (region  3).  This  phase  is  shown  in  plots 
c  and  d  of  figures  43  and  44.  Additional  tensile  failure  in  the  3-direction  occurs  on  the  face 
opposing  the  fixed  portion,  shown  in  Region  4  of  plots  c  and  d.  The  final  phase  is  compressive 
1 -direction  failure  in  the  0°  plies  starting  in  Region  1  of  plot  c,  expanding  to  Region  3  in  plots  d 
and  e,  and  also  occurring  in  Region  4  in  plots  d  and  e. 

Figures  36-44  display  the  information  that  is  needed  to  gain  insight  into  the  ply-level  failure 
behavior  of  the  sample.  This  insight  is  used  to  make  changes  to  the  layup  of  the  sample  to  meet 
the  objectives  of  this  case  study.  The  current  [0790°]  s  layup  does  not  satisfy  both  of  the 
objectives  of  the  case  study.  The  displacement  at  point  A  is  0.0585  m,  which  exceeds  the  target 
value  of  0.05  m.  The  areas  of  global  XY  shear  failure  (figure  39)  and  global  X-direction  failure 
(figure  41)  are  large  and  limiting  these  areas  is  the  other  important  objective. 

4.4.3  Design  Iteration  No.  1 

The  original  [0°/90°]s  layup  of  the  laminate  is  insufficient  in  both  the  displacement  objective  and 
limiting  the  extent  of  stiffness  loss.  The  design  is  first  changed  to  a  [02/±45°]s  layup.  The  0° 
plies  of  the  laminate  are  essential  in  providing  enough  longitudinal  stiffness  to  satisfy  the 
displacement  objective.  Using  ±45°  plies  is  aimed  at  addressing  the  area  of  shear  failure  along 
the  right  fixed  boundary  condition.  The  fiber  direction  of  these  plies  is  aligned  more  in  the 
global  X-direction  than  the  90°  plies,  adding  to  the  stiffness  of  the  structure  in  this  direction.  In 
the  original  design,  the  90°  plies  only  provided  minimal  stiffness  in  the  X-direction. 

In  order  to  determine  if  the  new  layup  has  improved  upon  the  original  design,  the  stiffness  ratios 
are  first  examined.  In  figure  46,  the  MINCIJ  of  the  new  design  is  compared  to  the  original.  The 
new  design  does  not  have  areas  with  a  complete  loss  in  stiffness  in  any  direction  (dark  blue 
areas).  The  new  design  has  only  a  small  area  in  the  lower  right  corner  where  a  minimum 
stiffness  ratio  has  not  changed  (red  area)  from  its  original  value.  In  the  original  design,  the  area 
of  no  loss  of  original  stiffness  is  large. 

The  majority  of  the  sample  in  right  plot  of  figure  46  has  some  stiffness  loss.  The  plot  shown  in 
figure  47  is  used  to  determine  in  which  direction  these  losses  are  the  greatest.  In  this  contour 
plot,  the  areas  in  red  signify  the  Y -direction  (C YYR)  and  the  areas  in  green  signify  the  X- 
direction  (CXXR),  the  blue  area  indicates  that  the  minimum  stiffness  ratio  was  not  below  80%. 
While  the  stiffness  in  the  X-direction  is  critical  to  the  response  of  the  structure,  stiffness  loss  in 
the  Y-direction  is  not  as  critical  because  there  is  no  loading  in  this  direction  though  there  is  a 
Poisson’s  effect. 
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Figure  46.  Comparison  of  minimum  stiffness  ratio  MINCIJ  [0°/90°]s  (left)  vs.  [02/±45°]s  (right). 
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Figure  47.  Plot  of  the  direction  of  minimum  Cy  ratio  MINCIJN  for 
[02/±45]s  layup. 


A  comparison  of  the  individual  stiffness  ratios  of  Cxx,  Cyy,  Czz,  and  Cxy  will  determine  the 
benefits,  if  any,  of  switching  to  this  layup.  The  comparison  plot  of  Cxx  in  figure  48  shows  that 
while  the  stiffness  loss  in  the  X-direction  of  the  [02/±45]s  sample  covers  a  larger  area  than  the 
original  design,  the  magnitude  of  stiffness  loss  in  these  areas  is  less.  In  this  plot,  the  areas  in 
light  blue  have  a  higher  stiffness  ratio  than  those  in  dark  blue. 
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Figure  48.  Plot  of  Cxx  stiffness  ratio  (CXXR)  comparing  [0°/90°]s  (left)  to  [02/±45]s  (right). 


Comparing  the  Cyy  stiffness  ratios  of  the  sample  in  figure  49,  it  is  seen  that  the  new  layup  has  a 
large  area  where  stiffness  has  been  reduced.  This  reduced  stiffness  is  a  result  of  the  0°  plies 
failing  in  the  transverse  matrix  direction,  which  causes  a  small  reduction  in  stiffness  in  this  area. 
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Figure  49.  Plot  of  Cyy  stiffness  ratio  (CYYR)  comparing  [0°/90°]s  (left)  to  [02/±45]s  (right). 


Similarly,  the  comparison  of  Czz  ratio  shown  in  figure  50  illustrates  that  the  complete  failure  of 
some  elements  in  the  Z-direction  has  been  replaced  by  slight  softening  over  the  entire  structure. 
The  response  of  the  sample  in  this  direction  is  determined  by  the  global  XZ  Poisson  ratio  and  the 
new  layup  has  less  coupling  in  this  direction. 
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Figure  50.  Plot  of  Czz  stiffness  ratio  (CZZR)  comparing  [0°/90°]s  (left)  to  [02/±45]s  (right). 


Finally,  an  examination  of  the  global  XY  shear  stiffness  ratio  between  the  original  and  modified 
layups  is  shown  in  figure  5 1 .  The  modified  layup  shows  dramatic  improvements  in  limiting  the 
extent  and  severity  of  shear  failure  along  the  fixed  boundary. 

The  results  of  this  first  design  iteration  are  a  reduction  in  the  area  and  magnitude  of  XY  shear 
stiffness  loss,  a  reduction  in  magnitude  of  Z  stiffness  loss,  more  widespread  but  less  severe  Y 
stiffness  loss,  and  a  decrease  in  the  magnitude  of  stiffness  loss  in  the  X-direction.  The 
displacement  of  point  A  for  this  layup  is,  however,  0.0647  m,  which  is  worse  than  the  first  layup 
tested.  While  there  is  improved  performance  in  the  stiffness  objectives,  the  displacement 
objective  is  still  not  met  for  this  layup. 
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Figure  51.  Plot  of  CXY  stiffness  ratio  (CXYR)  comparing  [0Y90°]S  (left)  to  [02/±45]s  (right). 


4.4.4  Design  Iteration  No.  2 

For  the  next  design  iteration,  the  stiffness  in  the  longitudinal  direction  of  the  structure  must  be 
improved  to  meet  the  displacement  objective.  A  laminate  in  a  [02/+35°]s  layup  is  chosen  because 
the  ±35°  plies  are  needed  to  address  the  XY  shear  failures  that  occur  along  the  fixed  boundary  of 
the  part.  These  plies  also  increase  the  global  X-direction  stiffness  of  the  sample  since  their  fibers 
are  oriented  more  in  this  direction  than  the  previous  design. 

The  first  condition  to  check  is  the  displacement  objective.  Examining  the  results  of  this  sample 
shows  that  the  X  displacement  at  point  A  is  once  again  not  within  the  desired  objective,  at 
0.0585  m.  An  investigation  into  the  stiffness  ratio  parameters  shows  that  there  is  an  improved 
performance  in  the  X  stiffness,  slightly  less  desirable  performance  in  the  Y  and  Z  directions,  and 
similar  performance  in  the  XY  shear  stiffness  when  compared  to  the  first  design  iteration. 
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4.4.5  Design  Iteration  No.  3 

The  next  layup  needs  to  increase  the  stiffness  of  the  laminate  in  the  X-direction  while  continuing 
to  prevent  shear  stiffness  loss.  A  laminate  in  a  [02/±25°]s  layup  satisfies  both  these  criteria  and  is 
used  for  this  design  iteration. 

The  first  check  for  the  layup  is  to  compare  the  displacement  at  point  A  to  the  targeted  value.  A 
displacement  of  0.0475  m,  less  than  the  targeted  value,  is  recorded  at  this  point.  With  this 
objective  satisfied,  an  examination  of  the  four  key  stiffness  ratios  is  needed  to  confirm  that  the 
performance  of  the  layup  meets  the  other  objective. 


Like  the  first  design  iteration,  an  examination  of  the  MINCIJ  of  this  layup  is  compared  to  the 
original  (figure  52).  The  majority  of  the  sample  with  the  new  layup  has  a  significant  loss  in 
stiffness  in  at  least  one  direction.  However,  in  contrast  with  the  original  layup,  this  layup  does 
not  have  large  areas  of  complete  stiffness  loss.  The  new  layup  has  only  one  element,  near  the  top 
of  the  fixed  edge,  with  complete  stiffness  loss. 


SDV15 
(Avg:  75%) 


+1. 000e+00 
+9. 583e-01 
+9.167e-01 
+8. 750e-01 
+8. 333e-01 
+7. 917e-01 
+7. 500e-01 
+7. 083e-01 
+6. 667e-01 
+6. 250e-01 
+5. 833e-01 
+5. 417e-01 
+5. 000e-01 
+4. 583e-01 
+4. 167e-01 
+3. 750e-01 
+3.333e-01 
+2. 917e-01 
+2.500e-01 
+2.083e-01 


i 


+1. 667e-01 
+1.250e-01 
+8. 333e-02 
+4.167e-02 
+0. 000e+00 


Figure  52.  Comparison  of  minimum  stiffness  ratio  MINCIJ,  [0°/90°]s  (left)  vs.  [02/±25°]s  (right). 
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The  plot  shown  in  figure  53  displays  that  the  Y-direction,  here  shown  in  pale  green,  is  the 
predominate  loss  of  stiffness  direction.  Stiffness  loss  in  the  Y-direction  is  not  critical  because 
there  is  no  loading  in  this  direction. 
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Figure  53.  Plot  of  the  direction  of  minimum  Qj  ratio  MINCIJN  for  [02/±25]s  layup. 
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A  comparison  of  the  individual  stiffness  ratios  of  Cxx,  Cyy,  C/z,  and  CXy  is  once  again 
performed.  The  comparison  plot  of  Cxxin  figure  54  shows  that  the  stiffness  loss  in  the  X- 
direction  has  been  minimized.  In  this  plot,  the  areas  in  blue  have  a  high  stiffness  loss  and  the 
areas  in  red  have  no  stiffness  loss. 
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Figure  54.  Plot  of  Cxx  stiffness  ratio  (CXXR)  comparing  [0°/90°]s  (left)  to  [02/±25]s  (right). 


Comparing  the  Cyy  stiffness  of  the  sample  in  figure  55,  it  is  seen  that  the  new  layup  has  a  large 
area  where  stiffness  has  been  reduced.  Because  stiffness  in  this  direction  is  not  critical  to  the 
overall  response  of  the  sample,  this  result  is  not  a  concern. 
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Figure  55.  Plot  of  CYy  stiffness  ratio  (CYYR)  comparing  [0°/90°]s  (left)  to  [02/+25]s  (right). 


Similar  to  the  first  design  iteration,  the  comparison  of  Czz  shown  in  figure  56  illustrates  that  the 
complete  failure  of  some  elements  in  the  Z-direction  has  been  replaced  by  slight  softening  over 
the  entire  structure. 
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Figure  56.  Plot  of  Czz  stiffness  ratio  (CZZR)  comparing  [0°/90°]s  (left)  to  [02/±25]s  (right). 


Finally,  an  examination  of  the  XY  shear  stiffness  ratio  between  the  original  and  modified  layups 
is  shown  in  figure  57.  The  modified  layup  shows  improvements  in  limiting  the  extent  and 
severity  of  shear  failure  along  the  fixed  boundary  when  compared  to  the  original  layup.  This 
layup  does  not  perform  as  well  as  the  first  design  iteration,  but  it  does  eliminate  the  shear  failure 
region. 

This  layup  is  able  to  satisfy  both  objectives  set  at  the  start  of  the  case  study.  With  a  displacement 
of  0.0475  at  point  A,  the  sample  is  within  the  targeted  value  of  0.05.  The  modifications  to  the 
layup  limit  the  areas  of  total  stiffness  loss  of  the  structure  in  the  global  X  and  XY  shear 
directions. 
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Figure  57.  Plot  of  CXY  stiffness  ratio  (CXYR)  comparing  [0°/90°]s  (left)  to  [02/±25]s  (right). 

This  case  study  demonstrates  the  utility  of  the  stiffness  ratio,  last  failed  mode,  and  last  failed  ply 
LAMPATNL  output  parameters.  These  parameters  helped  to  improve  the  understanding  of  the 
complex  interactions  between  progressive  ply  failure  and  global  stiffness  reduction  in  the 
laminate.  Areas  of  critical  stiffness  loss  were  easily  identified  using  the  linear  scale  of  the 
minimum  stiffness  ratio  parameter.  The  progression  of  ply  failures,  determined  from  the  design 
methodology,  provided  a  clear  understanding  of  the  response  of  the  structure.  Modifications  to 
the  [0°/90°]s  layup  were  made  based  on  the  insight  gained  from  the  mode  and  ply  failure  history. 
These  modifications  enabled  the  designer  to  quickly  achieve  the  performance  objectives  of  the 
analysis  using  only  three  design  iterations. 

4.5  Case  Study  No.  2 

4.5.1  Description 

The  second  case  study  is  3-D  model  of  an  open  hole  sample  that  is  subjected  to  multiaxial,  in¬ 
plane  loading.  The  open  hole  configuration  was  selected  because  it  is  a  standard  example  used 
to  examine  progressive  failure  models  for  composites  {18,  19,  21-25,  27,  29,  30).  This  problem 
is  also  common  in  composite  applications  and  the  effects  on  the  structural  response  due  to 
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changes  in  material,  layup,  and  width-to-diameter  ratio  are  important  to  understand.  A  diagram 
of  the  sample  is  shown  in  figure  58.  The  diameter  of  the  hole  is  1  and  the  width  of  the  sample  in 
both  the  X-  and  Y-directions  is  5,  resulting  in  width-to-diameter  ratio  of  5. 


Figure  58.  Diagram  of  open  hole  sample. 

The  original  laminate  is  an  S-Glass/Epoxy  material  in  a  [0°/90o/±45°]s,  or  quasi-isotropic,  layup. 
The  Ramberg-Osgood  parameters  and  maximum  strain  allowables  for  this  material  are  shown  in 
tables  6  and  7,  respectively  (30,  35).  The  mesh  for  this  model  contains  6768  elements  and  8905 
nodes.  The  FE  model  has  symmetry  boundary  conditions  along  the  left  face  and  bottom  face.  A 
uniform  pressure  load  of  300  MPa  is  applied  to  the  right  exterior  faces  in  the  positive  X 
direction,  creating  tensile  stresses  in  that  direction.  Another  pressure  load  of  150  MPa  is  applied 
to  the  top  exterior  faces  in  the  negative  Y  direction,  creating  compressive  stresses  in  that 
direction. 

The  objectives  in  this  design  case  study  are  to  limit  the  displacement  of  the  point  at  the  upper- 
right  comer  of  the  structure  to  less  than  0.05  m  in  the  positive  X-direction  and  0.04  m  in  the 
negative  Y-direction.  Based  upon  the  performance  of  the  quasi-isotropic  layup,  modifications  to 
the  laminate  may  be  needed  to  achieve  these  objectives.  Any  change  to  the  laminate  must 
maintain  the  same  material,  weight,  and  thickness  from  the  original  layup.  These  changes  are 
made  based  on  information  gained  from  the  process  outlined  in  the  design  methodology. 
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4.5.2  Results  of  Original  Design 

Following  the  procedure  of  the  design  methodology,  the  first  contour  plot  examined  is  a  plot  of 
the  minimum  Cy  ratio  MINCIJ,  shown  in  figure  59. 


Figure  59.  Contour  plot  of  minimum  Qj  stiffness  ratio  MINCIJ. 


The  yellow  and  bright  green  areas  of  the  sample  represent  a  stiffness  ratio  of  around  60%-70%, 
the  pale  blue  area  has  a  stiffness  ratio  of  10%-20%,  the  dark  blue  area  has  0%  stiffness  ratio,  and 
the  red  area  (2  elements)  is  near  100%  stiffness  ratio. 

The  multiaxial  loading  that  is  applied  to  this  sample  causes  stress  concentrations  at  both  sides  of 
the  open  hole  and,  as  expected,  the  stiffness  ratio  is  lowest  in  these  areas.  Almost  the  entire 
sample  has  a  30%  reduction  of  stiffness  in  at  least  one  direction.  To  determine  which  directions 
are  represented  in  figure  59,  a  plot  of  the  minimum  stiffness  ratio  number  MINCIJN  is  needed. 
This  plot  is  shown  in  figure  60. 
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Figure  60.  Contour  plot  of  the  direction  of  minimum  Cy  ratio  MINCIJN. 

In  this  figure,  light  blue  (MINCIJN=1)  represents  the  X-direction,  pale  green  (MINCUN=2) 
represents  the  Y-direction,  and  bright  green  (MINCUN=3)  represents  the  Z-direction.  The 
important  information  from  this  plot  is  that  the  area  of  stiffness  ratio  close  to  10%  at  the  top  of 
the  hole  corresponds  to  the  global  X-direction  and  the  area  of  complete  loss  in  stiffness  to  the 
right  of  the  hole  corresponds  to  the  Z-direction.  The  rest  of  the  sample  that  has  a  30%  reduction 
in  stiffness  is  a  combination  of  the  X-  and  Y-directions. 

The  reduction  of  stiffness  at  the  top  and  right  of  the  hole  are  an  indication  that  failure  probably 
has  occurred  in  these  areas.  The  contour  plot  of  the  number  of  failures  in  figure  61  shows  that 
these  areas  have  the  highest  amount  of  ply  failures  in  the  sample.  As  a  result,  an  investigation 
into  all  of  the  critical  stiffness  ratios  is  needed  because  the  stiffness  reduction  may  have  occurred 
in  multiple  directions  in  addition  to  what  is  indicated  in  figures  59  and  60. 
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Figure  61.  Contour  plot  of  NOF,  case  study  no.  2. 


The  multiaxial  loading  combined  with  XZ  and  YZ  Poisson  ratio  effects  contributes  to  strains  in 
the  Z-direction  of  the  sample.  The  softening  and  failure  of  the  laminate  in  the  through-thickness 
direction  is  shown  in  figure  62.  The  majority  of  this  sample  has  a  15%  reduction  in  stiffness,  the 
area  above  the  open  hole  has  been  reduced  by  35%,  and  the  area  to  the  right  of  the  hole  has  a 
complete  loss  in  stiffness.  This  stiffness  is  not  important  to  the  overall  response  of  the  sample, 
but  its  contribution  to  the  minimum  stiffness  ratio  shown  in  figures  59  and  60  could  mask  the 
reduction  in  other  important  stiffness  ratios. 
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Figure  62.  Plot  of  Czz  stiffness  ratio  (CZZR),  case  study  no.  2. 


The  next  stiffness  ratio  examined  is  the  Cxx  ratio  that  is  shown  in  figure  63.  The  three  areas  to 
note  in  this  plot  are  the  top  of  the  hole,  the  right  of  the  hole,  and  the  remainder  of  the  sample.  At 
the  top  of  the  hole,  the  Cxx  stiffness  has  been  reduced  to  15%  of  its  original  value.  To  the  right 
of  the  hole  and  in  the  rest  of  the  sample,  this  stiffness  ranges  from  50%  to  75%  of  its  original 
value.  It  is  important  to  note  that  in  the  area  to  the  right  of  the  hole,  this  stiffness  ratio  has  only 
been  reduced  to  50%  of  its  original  stiffness.  This  suggests  that,  although  some  plies  in  this 
region  may  have  failed  in  the  global  X-direction,  the  0°  plies  in  this  direction  have  not  failed  in 
fiber  tension.  The  ply  failure  history  of  this  sample  is  documented  later  in  this  section. 
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Figure  63.  Plot  of  Cxx  stiffness  ratio  (CXXR),  case  study  no.  2. 


The  plot  of  the  Cyy  stiffness  ratio  in  figure  64  has  similar  characteristics  to  the  previous  plot.  At 
the  right  of  the  hole,  the  Cyy  stiffness  has  been  reduced  to  15%  of  its  original  value.  This 
significant  drop  in  stiffness  is  not  apparent  when  looking  at  the  minimum  stiffness  ratio  plots. 

To  the  top  of  the  hole  and  in  the  remainder  of  the  sample,  this  stiffness  ranges  between  50%  and 
75%  of  its  original  value. 
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Figure  64.  Plot  of  Cyy  stiffness  ratio  (CYYR),  case  study  no.  2. 

The  final  stiffness  ratio  that  is  important  to  the  response  of  the  structure  is  the  Cxy  stiffness  ratio, 
which  represents  the  XY  shear  stiffness.  In  figure  65,  the  majority  of  the  sample  has  less  than 
5%  stiffness  loss  in  this  direction.  Areas  in  close  proximity  to  the  hole  in  the  sample  have  been 
reduced  to  65%  of  the  original  shear  stiffness. 
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Figure  65.  Plot  of  Cxy  stiffness  ratio  (CXYR),  case  study  no.  2. 


The  progression  of  failure  through  the  sample  must  be  understood  before  making  changes  to  the 
layup.  A  progression  plot  of  the  last  failed  mode  is  shown  in  figure  66.  For  this  plot,  dark  blue 
indicates  no  failure  (FMODE=0),  the  blue  areas  (FMODE=l)  indicate  1 -direction  tensile  failure, 
the  light  blue  areas  (FMODE=2)  indicate  1 -direction  compressive  failure,  the  areas  in  pale  green 
(FMODE=3)  are  2-direction  tensile  failure,  bright  green  areas  (FMODE=5)  have  tensile  failure 
in  the  3-direction,  and  red  (FMODE=9)  indicates  a  failure  mode  of  12  shear. 

A  progression  plot  of  the  last  failed  ply  is  shown  in  figure  67.  In  this  figure,  grey  areas  indicate 
no  failure  recorded  (FPLY=0),  yellow  areas  refer  to  the  90°  plies  (FPLY=7),  red  areas  to  the  0° 
plies  (FPLY=8),  blue  areas  to  the  -45°  plies  (FPLY=5),  and  finally  pale  green  areas  refer  to  +45° 
plies  (FPLY=6). 
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Figure  66.  Plot  of  progression  of  last  failed  mode  FMODE,  case  study  no.  2. 
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Figure  67.  Plot  of  progression  of  last  failed  ply  FPLY,  case  study  no.  2. 


The  quasi-isotropic  layup  of  the  laminate  makes  the  interpretation  of  these  plots  more  difficult 
than  the  first  case  study.  In  order  to  better  describe  the  progression  of  failure  in  the  sample,  five 
regions  are  identified  and  highlighted  on  the  sample  in  figure  68.  The  regions,  shown  in 
figure  68  on  a  plot  of  FMODE  at  the  end  of  the  analysis,  are  used  to  refer  to  areas  in  the 
progression  plots  of  figure  66a-g  and  figure  67a-g. 


Region  5 


Region  2 


Region  1 


Region  4 


Region  3 


Figure  68.  Diagram  of  regions,  case  study  no.  2. 


Failure  begins  in  this  sample  in  region  1,  the  top  portion  of  the  hole,  with  transverse  (2-direction) 
tensile  failure  in  the  90°  plies  (see  plot  a  of  figures  66  and  67)  and  propagates  across  almost  the 
entire  sample,  shown  in  the  yellow  areas  of  figure  67  plots  b  and  c.  The  next  failures  again  start 
in  region  1  during  plot  b  and  are  transverse  tensile  failure  in  the  ±45°  plies.  These  failures 
expand  outward  to  region  2  throughout  the  rest  of  the  loading  as  shown  in  the  green  and  blue 
areas  in  regions  1  and  2  of  figure  67  plots  c  to  g.  An  area  of  tensile  3-direction  failure  in  all  plies 
starts  in  region  3,  to  the  right  of  the  hole,  in  plot  c  and  expands  to  encompass  this  region  from 
plots  d  to  g.  Also  in  this  region,  both  compressive  transverse  failure  in  the  0°  plies  and  tensile 
transverse  failure  in  the  90°  plies  occur  from  plots  d  to  g.  A  pocket  of  transverse  tensile  failure 
in  the  ±45°  plies  (blue  area  in  region  4  of  figure  67e),  followed  by  compressive  transverse  failure 
in  the  0°  plies  (red  area  in  region  4  of  figure  67f),  arises  and  expands  in  region  4.  These  areas  of 
failure  expand  into  region  5  through  plots  f  and  g.  Finally,  12  shear  failure  occurs  in  the  ±45° 
plies,  which  is  the  red  area  in  region  4  of  figure  66g. 
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The  quasi-isotropic  layup  results  in  a  displacement  of  0.074  m  in  the  positive  X-direction  and 
0.046  m  in  the  negative  Y-direction.  These  values  are  not  within  the  desired  displacement 
objective  in  both  directions,  so  changes  to  the  layup  must  be  made  to  address  these  problems.  In 
addition,  there  is  a  large  area  of  failure  in  this  sample  and  minimizing  this  area  is  essential  in 
meeting  the  displacement  objective. 

4.5.3  Design  Iteration  No.  1 

Any  changes  to  the  original  layup  must  simultaneously  improve  stiffness  in  the  global  X-  and 
Y-directions,  because  each  displacement  objective  has  not  been  met.  To  achieve  this,  no 
additional  plies  can  be  added  to  the  new  layup.  In  the  [0°/90°/±45°]s  layup,  the  0°  plies  are 
essential  to  the  response  in  the  X-direction  as  are  the  90°  plies  for  the  response  to  the  Y-direction. 
Changing  the  remaining  ±45°  plies  to  additional  0°  and  90°  plies  would  help  achieve  the 
displacement  objectives,  but  there  is  XY  shear  failure  already  in  the  sample  and  losing  the  cross- 
plies  would  increase  this  failure.  A  change  to  the  ±45°  is  needed  to  address  either  the  X-direction 
or  the  Y-direction.  Since  the  displacement  in  the  X-direction  is  farther  from  its  objective, 
orienting  the  fibers  more  in  this  direction  will  be  beneficial. 

A  [0°/90°/±30°]s  layup  is  tested  in  this  design  iteration.  A  check  of  the  displacements  at  the  top 
right  comer  of  the  sample  shows  improvements  in  both  directions.  The  displacement  in  the 
positive  X-direction  is  0.052,  greater  than  the  targeted  value.  The  displacement  in  the  negative 
Y-direction  is  0.0405,  slightly  greater  than  the  targeted  value. 

This  analysis  shows  a  decrease  in  the  magnitude  of  the  displacement  in  the  Y-direction.  This 
result  is  unexpected  because  cross-ply  fibers  are  oriented  less  in  this  direction  when  compared  to 
the  original  layup,  yet  the  displacement  in  this  direction  has  decreased.  In  the  original  analysis, 
there  were  several  areas  of  the  sample  that  had  transverse  tensile  failure  in  the  ±45°  plies.  These 
failures  were  caused  by  the  tensile  loading  in  the  X-direction  of  the  sample  and  result  in  the 
2-direction  stiffness  in  these  plies  to  be  discounted.  The  overall  stiffness  of  the  structure  in  the 
Y-direction  is  softened,  causing  more  deflection. 

A  comparison  between  the  original  layup  and  the  current  layup  of  the  Cyy  stiffness  ratios  at  the 
end  of  loading  is  shown  in  figure  69.  The  initial  stiffness  of  the  new  layup  in  the  Y-direction  is 
less  than  that  of  the  original  layup,  but  the  stiffness  ratio  at  the  end  of  loading  is  greater.  The 
area  at  the  top  of  the  open  hole  retains  30%  more  stiffness  when  compared  to  original  layup. 

The  bulk  of  the  sample  has  a  stiffness  ratio  ranging  from  70%  to  95%,  compared  to  50%  to  75% 
in  the  original  layup.  The  area  to  the  right  of  the  hole  (region  3)  exhibits  more  stiffness  loss  than 
the  original  layup,  a  5%  ratio  compared  to  a  15%  ratio. 
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Figure  69.  Plot  of  CYY  stiffness  ratio  (CYYR)  comparing  [0°/90°/+45°]s  (left)  to  [0°/90o/±30°]s  (right). 

The  prevention  of  stiffness  loss  in  the  new  layup  comes  from  limiting  the  transverse  tensile 
failure  in  the  ±30°  plies,  illustrated  below  in  figure  70.  The  plot  shows  the  last  failed  ply  at  the 
end  of  loading.  The  areas  in  dark  blue  and  pale  green  represent  cross-plies  that  have  failed.  The 
corresponding  failure  mode  in  these  areas  is  transverse  (2-direction)  tensile  failure.  In  the  left 
plot  of  figure  70,  the  transverse  cross-ply  failure  region  is  large,  expanding  across  the  lower 
portion  of  the  sample  (region  5)  and  up  from  the  top  of  the  hole  (regions  1  and  2).  In  the  right 
plot  of  figure  70,  the  area  of  cross-ply  failure  is  only  a  small  portion  on  the  top  edge  of  the  hole 
(region  1),  significantly  smaller  than  the  original  layup. 


Figure  70.  Plot  of  last  failed  ply  (FPLY)  comparing  [0°/907+45°]s  (left)  to  [0°/90°/+30°]s  (right). 


90 


4.5.4  Design  Iteration  No.  2 

Building  on  the  results  and  insights  of  the  first  design  iteration,  additional  changes  need  to  be 
made  to  the  layup  to  meet  the  desired  performance  objectives.  The  displacement  of  the  comer 
point  in  the  negative  Y-direction  must  be  below  0.04  m,  so  improving  the  stiffness  and  reducing 
progressive  failure  in  this  direction  is  once  again  needed.  Based  on  the  first  iteration,  aligning 
the  fibers  further  into  the  X-direction  will  prevent  transverse  tensile  failure  in  these  plies, 
minimizing  stiffness  loss  in  the  Y-direction,  and  decreasing  the  magnitude  of  the  Y- 
displacement. 

The  layup  used  for  this  design  iteration  is  [0°/90o/±20°]s,  which  increases  the  fiber  orientation  of 
the  plus/minus  plies  by  10°  over  the  previous  iteration.  The  resulting  displacements  of  the  upper 
right  comer  of  the  sample  are  0.041  m  in  the  positive  X-direction  and  0.0342  m  in  the  negative 
Y-direction.  Both  displacements  are  within  their  respective  target  values,  so  this  layup  has 
achieved  the  performance  objectives.  Comparisons  of  the  critical  stiffness  ratios  will  display  the 
improvements  in  reducing  the  progressive  failure  within  the  structure. 

The  comparison  of  the  minimum  Qj  stiffness  ratio  MINCU  is  shown  in  figure  71.  When 
compared  to  the  original  layup,  the  improvements  are  not  apparent.  The  light  blue  area 
representing  a  15%  stiffness  ratio  in  the  new  layup  (figure  71  right)  is  smaller  and  confined  close 
the  edge  of  the  hole  (region  1).  The  area  of  complete  stiffness  loss  to  the  right  of  the  hole 
(region  3)  is  larger  in  the  new  layup.  Throughout  the  rest  of  the  sample,  the  minimum  stiffness 
ratios  are  greater  than  that  of  the  original  design. 


Figure  71.  Comparison  of  minimum  stiffness  ratio  (MINCU),  [07907+45°]s  (left)  to  [07907+20°]s  (right). 
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A  contour  plot  of  the  direction  of  minimum  Qj  stiffness  ratio  MINCIJN  is  shown  in  figure  72. 
The  areas  in  red  signify  the  XY  shear  ratio,  pale  green  signifies  the  Y-direction,  light  blue 
(1  element)  is  the  X-direction,  and  bright  green  is  the  Z-direction.  Dark  blue  areas  have  a 
minimum  stiffness  ratio  greater  than  80%  and  the  direction  is  not  recorded. 
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Figure  72.  Plot  of  the  direction  of  minimum  Qj  ratio  MINCIJN  for  [07907+20°] s  layup. 


The  plot  of  Cxx  stiffness  ratio  in  figure  73  shows  the  improved  performance  of  the  new  layup  in 
this  direction.  With  the  ±20°  fibers  aligned  more  in  the  X-direction  than  the  quasi-isotropic 
layup,  the  new  layup  is  expected  to  have  less  softening  and  failure  in  this  direction.  The  area  of 
85%  stiffness  loss  in  region  1,  the  top  of  the  hole,  for  the  original  layup  has  been  reduced  to  one 
element  along  the  edge  of  the  hole  in  the  new  layup.  The  rest  of  the  sample  has  become  just  5% 
to  25%  compliant,  an  improvement  on  the  25%  to  40%  compliance  of  the  original  layup. 
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Figure  73.  Plot  of  Cxx  stiffness  ratio  (CXXR)  comparing  [07907+45°]s  (left)  to  [07907+20°]s  (right). 


Like  the  first  design  iteration,  the  Cyy  stiffness  ratios  shown  in  figure  74  improve  upon  the 
original  layup.  Limiting  the  transverse  tensile  failure  in  the  ±20°  plies  improves  the  performance 
of  the  sample  in  the  Y-direction.  The  area  of  stiffness  loss  around  the  right  side  of  the  hole 
(region  3)  is  slightly  more  severe  than  the  original  layup,  but  the  improvement  in  stiffness  loss 
over  the  entire  sample  outweighs  this  drawback. 


Figure  74.  Plot  of  CYY  stiffness  ratio  (CYYR)  comparing  [0Y90Y+45°]S  (left)  to  [0Y90Y+20°]S  (right). 


The  Czz  stiffness  ratio  plot  in  figure  75  shows  improvement  from  the  original  layup  throughout 
the  entire  structure.  The  area  of  complete  stiffness  loss  in  right  plot  of  figure  71  is  shown  here  as 
stiffness  loss  in  the  Z-direction.  The  XZ  and  YZ  coupling  of  this  material  causes  strain  in  the 
Z-direction,  but  this  is  not  a  critical  direction  because  there  is  no  loading  in  this  direction. 
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Figure  75.  Plot  of  Czz  stiffness  ratio  (CZZR)  comparing  [07907+45 °]s  (left)  to  [07907+20°]s  (right). 

The  one  aspect  where  stiffness  loss  has  increased  from  the  original  design  to  the  new  design  is 
the  XY  shear  stiffness.  The  plot  of  the  Cxy  stiffness  ratio  shown  in  figure  76  illustrates  the 
extent  of  shear  stiffness  loss.  While  the  original  layup  had  very  little  area  of  concern,  the  new 
design  has  a  15%  stiffness  ratio  around  the  upper  edge  of  the  hole  (region  1)  and  a  stiffness  ratio 
between  50%  to  85%  throughout  the  rest  of  the  structure. 


Figure  76.  Plot  of  CXY  stiffness  ratio  (CXYR)  comparing  [07907+45°] s  (left)  to  [07907+20°]s  (right). 

The  [0°/90°/±20°]s  layup  achieves  the  displacement  performance  objectives  of  the  analysis. 
These  objectives  are  achieved  by  increasing  the  stiffness  of  the  structure  in  the  X-direction  and 
limiting  the  amount  of  transverse  tensile  failure  in  the  plus/minus  plies  of  the  layup. 
Comparisons  of  the  stiffness  ratios  between  the  original  and  new  layup  display  improvements  in 
reducing  the  amount  of  stiffness  loss  in  the  sample. 
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This  case  study  shows  that  the  approach  outlined  in  the  design  methodology  improves  the 
understanding  of  the  complex  material  nonlinearity  and  progressive  ply  failure  response  of  the 
open  hole  sample.  The  stiffness  ratio  parameters  provide  a  straightforward  measure  of  material 
degradation  and  reduce  the  complexity  of  the  results  when  compared  to  the  SF  output  parameter 
of  the  original  LAMPATNL  program.  Areas  of  critical  stiffness  loss  were  easily  identified  using 
the  minimum  stiffness  ratio  parameter.  The  failed  mode  and  ply  plots  for  this  sample  identified 
transverse  tensile  failure  in  the  ±45°  plies  of  the  quasi-isotropic  layup  as  a  cause  of  stiffness  loss 
in  the  global  Y-direction.  Modifications  to  the  [0°/90o/±45°]s  layup  were  made  to  address  this 
failure.  These  modifications  enabled  the  designer  to  quickly  achieve  the  performance  objectives 
of  the  analysis. 


5.  Conclusions 


The  development  of  a  methodology  to  design  and  analyze  thick- section  composite  structures 
using  nonlinear  ply  properties,  progressive  ply  failure,  and  a  homogenizing-dehomogenizing 
technique  (the  “smearing-unsmearing”  approach)  was  the  goal  of  this  research.  An  investigation 
into  the  ability  of  ABAQUS,  ANSYS,  and  MSC.Nastran  to  model  composite  materials  with 
anisotropic  nonlinear  behavior,  support  progressive  ply  failure  using  multiple  failure  criteria, 
homogenize  laminates  to  obtain  effective  properties,  and  dehomogenize  the  laminate  to  apply  ply 
based  failure  theories  was  conducted.  The  results  of  this  investigation  showed  that  none  of  the 
programs  contained  all  of  these  abilities.  Additional  composite  modeling  software,  such  as 
GENOA  and  Heilus:MCT,  were  also  researched  and  yielded  similar  conclusions. 

Presented  in  this  work  was  background  on  the  various  capabilities  and  limitations  of  the  “LAM” 
series  of  codes,  LAM3D,  LAM3DNLP,  LAMP  AT,  and  LAMPATNL.  These  codes  are  able  to 
model  thick-section  composite  structures  using  nonlinear  ply  properties,  progressive  ply  failure, 
and  the  “smearing-unsmearing”  approach.  In  the  “LAM”  codes,  ply-level  material  nonlinearity 
is  represented  by  using  the  Ramberg-Osgood  equation.  Progressive  ply  failure  uses  the  ply 
discount  method  in  combinations  with  ply-based  failure  criteria.  The  “smearing-unsmearing” 
approach  uses  the  theory  and  assumptions  of  Chou  (37)  to  calculate  laminate  and  ply-level 
stress-strain  response. 

Out  of  these  codes,  LAMPATNL  had  been  integrated  into  ABAQUS  using  a  UMAT  subroutine. 
LAMPATNL  is  the  LE  implementation  of  the  nonlinear,  analytical  model  LAM3DNLP.  Several 
modifications  were  made  to  the  LAM3DNLP  procedure  to  enable  it  to  function  in  the  LE 
environment  as  LAMPATNL.  The  modifications  include  changing  the  procedure  for  calculating 
ply-level  stresses  and  strains,  continuously  updating  of  the  Poisson’s  ratios  to  satisfy  material 
stability  conditions,  and  changing  the  post-failure  behavior  of  the  code.  Until  the  discussions  in 
this  research,  these  changes  were  never  documented  and  their  implications  were  never  discussed. 
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The  original  LAMPATNL  UMAT  was  not  able  to  accurately  track  the  failure  history  of  the 
elements  of  the  structure.  Also,  support  of  elements  was  limited  to  axisymmetric  and  plane 
strain  types.  A  major  contribution  from  this  research  was  adding  the  ability  to  track  the  failure 
history  by  element  number,  integration  point,  ply,  and  direction.  LAMPATNL  and  the  other 
“LAM”  codes  were  developed  to  analyze  thick-section  composite  structures,  so  3-D  analysis 
capabilities  were  needed.  Another  contribution  of  this  research  was  expanding  the  UMAT  to 
support  3-D  eight-node  linear  brick  elements. 

In  addition  to  the  undocumented  changes  to  LAM3DNLP  when  it  was  implemented  in 
LAMPATNL,  a  full  validation  of  the  stress-strain  response  output  of  this  UMAT  was  also  never 
provided.  With  four  materials  and  six  different  layups,  a  total  of  13  load  cases  were  simulated 
using  LAMPATNL  and  compared  to  results  generated  by  LAM3DNLP.  Both  linear  and 
nonlinear  ply  properties  were  tested  for  these  load  cases  to  ensure  full  compatibility  with  both 
LAM3D  and  LAM3DNLP.  None  of  the  example  cases  presented  for  validation  of  LAMPATNL 
deviate  from  the  values  produced  by  LAM3DNLP. 

5.1  Summary  of  Results 

LAMPATNL  had  the  original  output  parameters  of  safety  factor,  critical  mode,  and  critical  ply. 

It  was  extremely  difficult  to  determine  the  areas  and  directions  of  failure  in  the  structure  at  the 
end  of  analysis  using  these  parameters  alone.  An  example  was  provided  to  illustrate  the 
difficulty  of  using  the  safety  factor  to  analyze  a  composite  structure.  The  most  important 
contribution  of  this  research  was  the  stiffness  ratio  parameters  output  by  the  UMAT.  These 
parameters  improved  the  visualization  of  progressive  failure  and  stiffness  loss  and  significantly 
increased  the  utility  of  LAMPATNL. 

The  validated  LAMPATNL  UMAT  subroutine  was  then  used  to  analyze  thick-section  composite 
structures.  Another  unique  contribution  of  this  research  was  the  newly  developed  design 
methodology,  which  outlined  the  process  of  designing  a  structure  using  the  UMAT.  The 
methodology  used  stiffness  ratios  and  failure  progression  to  provide  insight  into  the  failure 
mechanics  of  the  laminate.  Changes  to  the  layup  of  the  composite  were  made  based  upon  the 
objectives  of  the  analysis  using  the  information  provided  by  the  design  methodology. 

Two  example  cases  were  used  to  evaluate  the  process  of  the  design  methodology.  The  first  case 
study  was  a  [0790°] s  layup  of  AS4/3501-6  subjected  to  compression  and  in -plane  shearing.  The 
objective  of  this  case  was  to  limit  the  displacement  in  the  X-direction.  This  objective  was 
achieved  by  increasing  the  orientation  of  fibers  in  X-direction  and  limiting  failure  in  the  sample. 
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The  resulting  design  for  this  case  was  a  [02/±25°]s  layup.  This  layup  was  able  to  achieve  the 
performance  objective  of  the  case  study  by  limiting  the  amounts  of  shear  and  tensile  failure  and 
increasing  the  longitudinal  stiffness  in  the  sample. 

The  second  case  study  was  an  open  hole  sample  of  S-Glass/Epoxy  in  a  quasi-isotropic  layup 
subjected  to  multiaxial  in-plane  loading.  The  failure  progression  of  this  structure  revealed  that 
transverse  tensile  failure  in  the  ±45°  plies  resulted  in  decreased  stiffness  in  the  Y-direction. 
Increasing  the  orientation  of  the  plus/minus  plies  into  the  X-direction  leads  to  improved  stiffness 
in  both  directions  due  to  the  limiting  of  transverse  tensile  failure  in  these  plies.  A  [0°/90o/±20°]s 
layup  was  determined  to  meet  both  objectives  of  the  analysis  while  keeping  the  ply  amount 
constant  and  therefore  not  adding  weight  to  the  structure. 

5.2  Future  Work 

The  documentation  and  validation  of  the  LAMPATNL  UMAT  subroutine  provides  a  basis  to 
start  analyzing  thick-section  composite  structures.  The  design  methodology  provides  a  template 
for  using  LAMPATNL  to  modify  nonlinear  composite  layups  using  progressive  ply  failure 
analysis.  The  two  case  studies  illustrated  the  application  of  UMAT  to  design  the  layup  of  two 
structures  under  complex  loading.  The  expansion  of  the  LAMPATNL  to  analyze  other  cases  of 
thick-section  composite  structures  is  one  subject  of  future  work.  In  addition  to  the  design  and 
analysis  of  these  structures,  experimental  validation  of  the  resulting  design  modifications  from 
LAMPATNL  is  also  an  important  step. 

Another  subject  of  future  work  from  this  research  is  expanding  the  LAMPATNL  user  material  to 
explicit  dynamic  loading.  Currently,  the  UMAT  developed  for  ABAQUS  is  only  compatible 
with  the  implicit  solver.  This  solver  is  used  to  simulate  the  response  of  structures  under  static  or 
implicit  dynamic  loading.  Modifications  to  the  LAMPATNL  code  that  interfaces  with 
ABAQUS  are  required  for  UMAT  to  function  in  the  explicit  solver  environment.  Lurther 
expansion  in  the  dynamic  modeling  field  necessitates  the  addition  of  rate  dependent  material 
behavior.  The  current  constitutive  model  used  for  LAMPATNL  is  strain  dependent.  Additional 
work  would  be  required  to  develop  and  validate  a  strain-rate-dependent  constitutive  model. 

The  ability  to  model  residual  thermal  stresses  that  resulted  from  processing  the  composite 
materials  was  present  in  the  LAM3D  code.  These  thermal  modeling  features  were  largely 
abandoned  in  the  transition  to  nonlinear  materials  and  were  similarly  absent  in  the  LAMP  AT  and 
LAMPATNL.  With  the  framework  of  this  capability  still  in  place,  thermal  modeling  features 
could  be  reintroduced  without  much  difficultly.  Strength  predictions  of  metal  matrix  composite 
structures  were  shown  to  be  effected  from  initial  thermal  stresses  in  the  structure  (40). 

Expanding  the  ability  of  UMAT  to  model  thermal  effects  is  a  goal  of  future  LAMPATNL 
development. 


97 


As  discussed  in  section  3.1,  there  are  changes  to  the  procedure  of  calculating  ply-level  stresses 
and  strains  from  LAM3DNLP  to  LAMPATNL.  In  that  section,  it  was  noted  that  these  changes 
do  not  affect  in-plane  stress-strain  predictions  but  may  have  implications  for  hybrid  composites 
and  laminates  exposed  to  out-of-plane  loading.  In  order  to  implement  LAM3DNLP  into  the  FE 
code,  it  was  required  to  not  calculate  the  out-of-plane  ply  strain  in  accordance  with  the  original 
theory.  Further  research  is  required  to  determine  and  discuss  the  effects,  if  any,  of  setting  all  ply 
strains  equal  to  total  laminate  strains. 
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